*FUT2* **Secretor Status Influences Susceptibility to VP4 Strain-Specific Rotavirus Infections in South African Children**

#### **Jaime MacDonald 1,2,\*, Michelle J. Groome <sup>3</sup> , Janet Mans <sup>2</sup> and Nicola Page 1,2**


Received: 31 July 2020; Accepted: 22 September 2020; Published: 27 September 2020

**Abstract:** Gastroenteritis is a preventable cause of morbidity and mortality worldwide. Rotavirus vaccination has significantly reduced the disease burden, but the sub-optimal vaccine efficacy observed in low-income regions needs improvement. Rotavirus VP4 'spike' proteins interact with FUT2-defined, human histo-blood group antigens on mucosal surfaces, potentially influencing strain circulation and the efficacy of P[8]-based rotavirus vaccines. Secretor status was investigated in 500 children <5 years-old hospitalised with diarrhoea, including 250 previously genotyped rotavirus-positive cases (P[8] = 124, P[4] = 86, and P[6] = 40), and 250 rotavirus-negative controls. Secretor status genotyping detected the globally prevalent G428A single nucleotide polymorphism (SNP) and was confirmed by Sanger sequencing in 10% of participants. The proportions of secretors in rotavirus-positive cases (74%) were significantly higher than in the rotavirus-negative controls (58%; *p* < 0.001). The rotavirus genotypes P[8] and P[4] were observed at significantly higher proportions in secretors (78%) than in non-secretors (22%), contrasting with P[6] genotypes with similar proportions amongst secretors (53%) and non-secretors (47%; *p* = 0.001). This suggests that rotavirus interacts with secretors and non-secretors in a VP4 strain-specific manner; thus, secretor status may partially influence rotavirus VP4 wild-type circulation and P[8] rotavirus vaccine efficacy. The study detected a mutation (rs1800025) ~50 bp downstream of the G428A SNP that would overestimate non-secretors in African populations when using the TaqMan®SNP Genotyping Assay.

**Keywords:** rotavirus; secretor status; histo-blood group antigens; VP4 genotypes; *FUT2*; susceptibility; vaccines

#### **1. Introduction**

Gastroenteritis is a preventable cause of morbidity and mortality worldwide, and the burden predominantly exists in high-risk populations such as children under the age of five years in low-income regions [1]. Rotavirus is the most frequent aetiology of diarrhoeal illness and death in children<5 years-old, and it was responsible for 29% of global diarrhoeal deaths occurring in this age group in 2016 [2].

The introduction of oral rotavirus vaccines in >100 countries worldwide has significantly reduced the burden of rotavirus diarrhoea and resulted in a 38% overall reduction in childhood diarrhoeal hospitalisations globally [3,4]. However, rotavirus vaccine efficacy appears to vary significantly between high-income (85–98%) and low-income (50–64%) countries [5]. Eliciting an adequate immune response to oral vaccines is multifactorial but may be limited in low-income settings due to impoverished living conditions and increased exposure to pathogens [3]. In addition, the passive transfer of rotavirus maternal

antibodies during breastfeeding can influence the immune response elicited by oral rotavirus vaccines in young children [4]. Understanding the factors that have contributed to an observed lower rotavirus vaccine efficacy in these settings may alleviate the burden of rotavirus-associated mortality in children.

Host genetic factors have recently been proposed to influence susceptibility to enteric pathogens. The excretion of soluble human histo-blood group antigen (HBGA) structures in gut mucosal surfaces determines a host's 'secretor status,' controlled by the human *FUT2* gene. Non-secretor phenotypes with an inability to express soluble HBGAs due to mutations in the *FUT2* gene (such as the prevalent G428A SNP; rs601338) are present globally in varying proportions. Higher proportions of non-secretor phenotypes are observed in African populations (~30%) than in Asian populations (~5%) [6,7].

Antigenic HBGA structures present in the body can act as receptors for various pathogens to bind during infection [8,9]. *FUT2* secretor status can modulate infection because it defines the presence (secretor) or absence (non-secretor) of HBGA attachment factors excreted in the gut. Susceptibility to enteric norovirus infection has been associated with secretor status, where non-secretor phenotypes have been found to display a natural resistance to GII.4 norovirus strains [10–12]. It has been proposed that variations in secretor status phenotypes and subsequent differences in host-defined susceptibility may contribute to the circulation of rotavirus strains in a similar mechanism [13].

Interactions between rotavirus particles and HBGA receptors present in the gut can occur via the VP4 (VP8\* subunit) 'spike' protein on the surface of the virion [14]. Evidence of rotavirus VP4 strain-specific binding patterns between HBGAs and prevalent strains (P[8], P[4], and P[6]) has recently been noted [14]. Rotavirus P-types have distinct VP4 morphology that determines the presence or absence of HBGA-binding interfaces, allowing for different mechanisms of binding and entry of rotavirus particles to occur [13]. Studies have shown that rotavirus genotypes P[8]- and P[4]-bound complex and soluble HBGAs abundant in secretors, as well as an increased susceptibility to infection with these rotavirus strains in secretors. Non-secretors with an absence of HBGAs in the gut have been found to display a natural resistance to P[8] and P[4] strains with VP4 HBGA-binding interfaces [15–17]. Variations in host-defined secretor status can therefore influence susceptibility to infection with different rotavirus strains.

Rotavirus P[8] genotypes are responsible for more than 80% of human wild-type infections globally [15]. However, rotavirus circulation in Africa differs in strain diversity and prevalence, with more frequent cases of P[6] strains, which have reached 26% of all rotavirus strains circulating in African populations [18]. The proportions of naturally resistant non-secretors may alter the circulation of rotavirus P-types compared to that in global populations.

The Rotarix® and RotaTeq® rotavirus vaccines both contain P[8]-based strains or reassortants, and they provide protection through the replication of live-attenuated vaccine strains in the gut to induce a local immune response [4]. Associations between host-defined secretor status and susceptibility to infection with specific rotavirus strains pose interesting questions surrounding the lowered efficacy of P[8]-based rotavirus vaccines observed in some regions [19,20]. Emerging research has alluded to this idea [21–23], including the influence of the related *FUT3* Lewis host genetic factor [24–27], but further investigations are required. These data have contributed to the evidence that host genetic factors such as secretor status can influence infections by pathogens including rotavirus, as well as that strain-specific interaction mechanisms may occur [14,15,28].

The aim of this study was to investigate *FUT2*-defined secretor status in South African children <5 years-old hospitalised with diarrhoea and to examine the association between a host's genetic secretor status and rotavirus-associated hospitalisations. Understanding the relationship between pathogens such as rotavirus and the genetics of a population may identify avenues for improvements in vaccine efficacy to reduce the burden of rotavirus gastroenteritis.

#### **2. Results**

Secretor genotypes were successfully determined for all 500 children selected for the study, and the total cohort comprised 65.8% (329) secretors with at least one functional *FUT2* allele and 34.2% (171) non-secretors with both *FUT2* alleles containing the G428A SNP.

Rotavirus-positive cases (RV+) comprised 74% (185/250) secretors (Se) and 26% (65/250) non-secretors, while rotavirus-negative controls (RV-) comprised 58% (144/250) secretors and 42% (106/250) non-secretors. The distributions of secretors versus non-secretors observed amongst cases and controls were significantly different (*p* < 0.001).

Information on rotavirus genotyping from the Rotavirus Sentinel Surveillance Program (RSSP) database [29,30] showed that the rotavirus-positive cases (n = 250) comprised 124 P[8] infections, 86 P[4] infections, and 40 P[6] infections (Supplementary Material). The proportions of secretors and non-secretors were compared amongst each VP4 strain within rotavirus-positive cases (Table 1). Rotavirus P[8] infections (79% secretors and 21% non-secretors) and P[4] infections (77% secretors and 23% non-secretors) had significantly different proportions of secretor phenotypes compared to P[6] infections (53% secretors and 47% non-secretors) (*p* = 0.001 and *p* = 0.006, respectively). When considered together, rotavirus P[8] and P[4] infections (78% secretors and 22% non-secretors) had significantly different proportions of secretor phenotypes compared to P[6] infections (53% secretors and 47% non-secretors) (*p* = 0.001).

**Table 1.** The distribution of secretors and non-secretors amongst VP4 genotypes P[8], P[4], and P[6] of rotavirus-positive cases (RV+; n = 250).


The Sanger sequencing of the exon 2 region of the *FUT2* gene conducted for 10% of the cohort confirmed the presence of either functional *FUT2* alleles or G428A SNP alleles for 91% (48/53) of analysed specimens. Sequences of the *FUT2* exon 2 region from 12 homozygous secretors (SeSe), 24 heterozygous secretors (Sese), and 17 homozygous non-secretors (sese) were obtained and compared to RT-PCR G428A genotyping results. Five discrepant results were observed in which heterozygous secretor (Sese) individuals (one functional *FUT2* allele and one allele containing the non-functional G428A SNP) genotyped by Sanger sequencing were incorrectly genotyped by RT-PCR as non-secretors (both alleles containing the G428A SNP). A commonality between these discrepant specimens was an SNP mutation (rs1800025) ~50 bp downstream of the G428A SNP (Figure 1).

**Figure 1.** Sequence alignment of five participants where Sanger sequencing and RT-PCR genotyping results were discrepant. (**a**) The G428A SNP location displaying all discrepant sequences containing the two peaks 'G' and 'A,' as represented by an 'R' annotation. (**b**) The mutation site (rs1800025) located ~50 base pairs downstream of the G428A SNP, common in all discrepant results.

#### **3. Discussion**

The results from this study indicate that secretors were more susceptible to rotavirus infection, and non-secretors seemed to display a natural resistance. The absence of HBGAs in the gastric mucosa of non-secretors appeared to reduce susceptibility to rotavirus, possibly by limiting the attachment stage of binding and entry during rotavirus infection [31]. Despite this observation, non-secretors were present amongst rotavirus-positive cases, indicating that HBGA attachment may not be the only mechanism for rotavirus binding and subsequent entry. Early studies on rotavirus binding and entry described sialic acid as an attachment factor for some animal strains [32]. Alternative binding receptors such as sialic acid or yet unknown mechanisms could explain the presence of rotavirus infection in non-secretor individuals in our study.

Studies have shown that rotavirus VP4 (VP8\*) binds to HBGAs in a strain-specific manner [13]. Xu and colleagues showed that P[8] and P[4] rotavirus strains similarly bound to complex HBGAs via a ββ binding domain, while more distantly related P[6] strains bound simple H-type 1 structures in a βα binding domain [14]. In our study, a higher proportion of secretors was observed in P[8] (78%) and P[4] (76%) rotavirus infections compared to P[6] infections (53%). This suggested that secretors were significantly more susceptible to P[8] and P[4] strains than to P[6] strains (*p* < 0.01), while non-secretors were more likely to be infected with P[6] strains. These strain-specific interactions may also influence the circulation of rotavirus strains within the South African population, as observed in other settings [15–17].

A correlation in the prevalence of rotavirus VP4 strains and HBGA genotypes suggested that the circulation of rotavirus may be partially modulated by their ability to bind to host-defined HBGA receptors. Globally, G1P[8] is the predominantly circulating rotavirus genotype, with ~74% of global strains containing the P[8] VP4 strain [18]. However, studies have shown that rotavirus strains in Africa are more diverse, with P[8] comprising 32% of rotavirus cases, P[4] comprising 13% of rotavirus cases, and P[6] comprising 26% of rotavirus cases [18]. In South Africa, P[6] strains were detected in 25% of rotavirus cases between 2003 and 2006, and they continue to circulate [30,33]. In this study, the higher proportion of non-secretors (34%), naturally resistant to P[8] and P[4] rotavirus infections, may explain the 16% detection of P[6] strains [17,34]. The *FUT2* genetics of a population may define the availability of host HBGA receptors for rotavirus infection, which could drive the epidemiology of rotavirus strain circulation in a region.

Discrepant results in Sanger sequencing revealed that five individuals were misclassified by RT-PCR as non-secretors (error rate 22.7%; 5/22), with sequencing identifying these five individuals as heterozygous secretors (Sese). The specimen sub-set comprised 58.5% secretors and 41.5% non-secretors based on RT-PCR genotyping, while the same specimens comprised 67.9% secretors and 32.1% non-secretors based on Sanger sequencing—an overall over-estimation of non-secretors of approximately 10%. This over-estimation of non-secretor genotypes is important to note for future studies, especially when using the TaqMan® SNP Genotyping Assay targeting the G428A SNP in an African population where non-secretors are frequent. The proportion of non-secretors (34%) observed in our cohort of 500 individuals correlated with other studies in African populations where higher frequencies of non-secretors were observed [35,36].

Misclassification by the commercial genotyping assay was hypothesised to be due to a mutation noted ~50 bp downstream of the G428A SNP position. The manufacturer confirmed that the mutation affected the primer binding of the reverse primer to the functional copy of the *FUT2* gene in the five heterozygous secretors, resulting in the absence of PCR product for the FAM-labelled probe (which detects the presence of the allele without the G428A SNP) to bind. Interestingly, the mutation was found in 9% of African populations compared to 2% in all populations in the 1000 genomes project [37]. Sanger sequencing remains an important tool to investigate host genetic factors such as secretor status, and further sequencing will be considered to examine the extent of the *FUT2* G514R mutation detected in this study.

Studies have indicated that secretor status can influence antibody titres to rotavirus [36], the incidence of gastrointestinal disease [38], and immune responses to rotavirus vaccines [28]. Rotarix® and RotaTeq® vaccines both contain P[8] vaccine constructs and require multiplication in intestinal cells to elicit local gut immunity [39,40]. The absence of HBGA attachment factors in non-secretors may reduce the replicative capacity of P[8] vaccine strains. The observation that non-secretors in Africa exhibit a natural resistance to wild-type P[8] strains may provide insights into the differences in vaccine efficacy across populations [5]. A study by Kazi and colleagues identified a link between the immune response to rotavirus P[8] vaccines and secretor status [28], and these associations have since been observed elsewhere [19,22,41]. Since patient sera were not collected as part of the RSSP, we could not investigate the direct effect of secretor status on rotavirus vaccine immune responses. Future studies investigating links between secretor status and variables such as vaccine immune responses, breastfeeding in young children, population genetics, and gut microbiome compositions, as well as alternative binding receptors for rotavirus entry, should be considered.

The limitations of this study include the small sample size of P[6] rotavirus cases available for further analysis (16%; 40/250). A larger sample size of rotavirus genotypes would be beneficial in confirming the relationship between specific rotavirus VP4 strains and secretor status. Another limitation of this study was the discordant results between RT-PCR genotyping and Sanger sequencing, resulting in the misclassification of heterozygous secretors by RT-PCR. Only 13% (22/171) of non-secretor genes were sequenced due to budget constraints, and additional funding will be sought to expand the sequencing of the *FUT2* gene of non-secretors in South Africa. A final limitation of this study was not including analysis of the related *FUT3* Lewis genes as it may also impact susceptibility to rotavirus infections. Future studies should consider the genetics of a cohort before utilising genotyping techniques, since alternative SNPs may be present which may skew results.

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

The South African RSSP enrolled children under the age of five years hospitalised for diarrhoea at various sites across South Africa (Protocol M091018, approved by the Human Research Ethics Committee (Medical) of the University of Witwatersrand). Diarrhoea was defined as three or more loose stools in past 24 h, with or without vomiting.

Informed consent was obtained from each child's parent or guardian prior to participation in the RSSP. Stool and dried blood spot (DBS) specimens were collected from enrolled participants, and each child's stool was screened as part of the RSSP for rotavirus group A (ProspecT™ Rotavirus Microplate Assay, Oxoid, Basingstoke, UK). Rotavirus-positive cases were genotyped using conventional RT-PCR methods and primers for G-specific and P-specific genotypes to determine the GxP[x] rotavirus strain [42].

This sub-study was conducted in accordance with the Declaration of Helsinki, and the project entitled "Investigation of secretor status, rotavirus VP4 genotypes, and gastrointestinal microbiomes in cases of diarrhoea in South Africa" (Protocol number 222/2018) was approved by the Research Ethics Committee, Faculty of Health Sciences, University of Pretoria, in May 2018.

For this study, children enrolled in the RSSP between 2009 and 2017 with available DBS specimens were identified, and rotavirus-negative cases (n = 250) were randomly selected. Rotavirus GxP[x] genotypes were previously determined as part of the RSSP [30], and the rotavirus-positive subset (n = 250) was selected to represent the major rotavirus VP4 genotypes (P[8], P[4], and P[6]), with cases and controls selected randomly where possible.

Secretor status was investigated using DBS specimens. DNA from DBS specimens was extracted using a QIAamp DNA Mini kit (Qiagen Inc., Valencia, CA, USA) according to the manufacturer's instructions with one modification prior to extraction. The manufacturer's protocol was modified to improve lysis by incubating DBS cards (~1 cm diameter) in a 200 µL buffer ATL overnight at 37 ◦C, instead of at 85 ◦C for 10 min. Following extraction, DNA was stored at −40 ◦C at the Centre for Enteric Diseases (Virology), National Institute for Communicable Diseases.

Secretor status was determined by detecting the presence or absence of the *FUT2* G428A SNP using a Predesigned TaqMan® SNP Genotyping assay (Life Technologies Corporation, CA, USA, supplied by Thermo Fisher Scientific, Carlsbad, CA, USA) in a 10 µL reaction volume according to the manufacturer's instructions [27,43].

The Sanger sequencing of 10% of the cohort *FUT2* genes was performed to ensure that alternative non-secretor-causing SNPs, which may be undetected by this assay, were absent. The specimens were selected to include all secretor genotypes, with a slight selection bias towards heterozygous secretors (n = 19) and non-secretors (n = 22) compared to homozygous secretors (n = 12), as well as a range of cycle threshold values (Ct range of 10–39) obtained during RT-PCR. The coding exon 2 region of the *FUT2* gene was amplified using the FUT2Ex2F and FUT2Ex2R primers [7], cleaned using an ExoSAP-IT™ PCR Product Cleanup protocol (Thermo Fisher), and sequenced using a BigDye™ Terminator v3.1 Cycle Sequencing kit (Applied Biosystems, Life Technologies, Waltham, MA, USA) on an Applied Biosystems 3500xL Genetic Analyzer instrument (Applied Biosystems). Sequences were aligned to a *FUT2* protein-coding reference sequence (NG\_007511.1:11987-13018 *Homo sapiens* fucosyltransferase 2 (*FUT2*), RefSeqGene on chromosome 19) (NCBI) using Molecular Evolutionary Genetics Analysis software version 7.0.26 (MEGA7).

The sequences of the FUT2 exon 2 region of 10% of the cohort were submitted to BankIt (National Center for Biotechnology Information, Bethesda, MD, USA), and the accession numbers are as follows: MW036696, MW036697, MW036698, MW036699, MW036700, MW036701, MW036702, MW036703, MW036704, MW036705, MW036706, MW036707, MW036708, MW036709, MW036710, MW036711, MW036712, MW036713, MW036714, MW036715, MW036716, MW036717, MW036718, MW036719, MW036720, MW036721, MW036722, MW036723, MW036724, MW036725, MW036726, MW036727, MW036728, MW036729, MW036730, MW036731, MW036732, MW036733, MW036734, MW036735, MW036736, MW036737, MW036738, MW036739, MW036740, MW036741, MW036742, MW036743, MW036744, MW036745, MW036746, MW036747, MW036748.

Statistical analyses using Chi-squared tests and univariate logistic regression models were performed using STATA version 14.0, where *p* < 0.05 was considered significant (StataCorp College Station, TX, USA).

#### **5. Conclusions**

Rotavirus susceptibility appeared to be influenced by secretor status in this study of South African children hospitalised with acute diarrhoea. Secretors expressing HBGAs in gut mucosal surfaces were more likely to be infected with rotavirus, specifically the P[8] and P[4] strains, compared to non-secretors. Non-secretors, with an absence of HBGAs in the gut, appeared to be less susceptible to rotavirus P[8] and P[4] infections compared to secretors—thus, the P[6] genotype was more frequent in these individuals. Interactions between rotavirus and secretor status could provide insights into the circulation of rotavirus strains amongst genetically diverse populations. Insights into the potential causes of altered rotavirus susceptibility and subsequent vaccine efficacy will aid in minimising the burden of disease. Diarrhoeal deaths are preventable, and secretor status may be an important host genetic factor to help understand and improve rotavirus disease prevention. Finally, the choice of assay for detecting or classifying secretor status in different populations should be carefully considered because the tools currently available all have pros and cons associated with their use.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/10/795/s1, File: Final DBS Cohort Results\_7.9.2020.

**Author Contributions:** Conceptualization, M.J.G., J.M. (Janet Mans) and N.P.; Data curation, J.M. (Jaime MacDonald); Formal analysis, M.J.G.; Funding acquisition, N.P.; Investigation, J.M. (Jaime MacDonald); Methodology, J.M. (Jaime MacDonald) and J.M. (Janet Mans); Project administration, N.P.; Resources, N.P.; Supervision, M.J.G., J.M. (Janet Mans) and N.P.; Validation, N.P.; Writing—original draft, J.M. (Jaime MacDonald); Writing—review & editing, M.J.G., J.M. (Janet Mans) and N.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Rotavirus Surveillance grant (GSK E-Track 200238) and a bursary from the Poliomyelitis Research Foundation (Grant 19/48).

**Acknowledgments:** I would like to acknowledge the Centre for Enteric Diseases team at the National Institute for Communicable Diseases for facilitating and supporting this research.

**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.

#### **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*

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

**Peter N. Mwangi <sup>1</sup> , Milton T. Mogotsi <sup>1</sup> , Sebotsana P. Rasebotsa <sup>1</sup> , Mapaseka L. Seheri <sup>2</sup> , M. Je**ff**rey Mphahlele <sup>3</sup> , Valantine N. Ndze <sup>4</sup> , Francis E. Dennis <sup>5</sup> , Khuzwayo C. Jere 6,7 and Martin M. Nyaga 1,\***


Received: 13 April 2020; Accepted: 18 May 2020; Published: 20 May 2020

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

**Keywords:** atypical strains; genome constellation; reassortment; rotavirus; whole-genome characterization

#### **1. Introduction**

Diarrhea persists as a leading infectious mortality cause in children under the age of five worldwide [1]. Group A rotavirus (RVA) is the primary viral etiologic agent for acute gastroenteritis in children under five years of age [2], resulting in annual mortality cases ranging from 122,322 to 215,757 with an estimated 81% reported in sub-Saharan Africa and Southeast Asia [3,4]. To combat RVA diarrhea, especially in countries with high RVA disease burden, the World Health Organization (WHO) recommends incorporation of RVA vaccines into the national immunization programs alongside other childhood vaccines [5]. The WHO has prequalified four vaccines (Rotarix®, GlaxoSmithKline, Rixenstart, Belgium; RotaTeq®, Merck & Co, USA; ROTAVAC®, Bharat Biotech, Hyderabad, India and ROTASIL®, Serum Institute of India, Pune, India) for global use [6]. Two vaccines (Rotavin-M1®, POLYVAC, Hanoi, Vietnam and Lanzhou lamb rotavirus, Lanzhou Institute of Biological Products, Lanzhou, China) have been approved for national use in Vietnam and China, respectively [7,8]. Human neonatal RVA vaccine (RV3-BB) and bovine human reassortant RVA vaccine candidates as well as neonatal and non-replicating injectable vaccines are in the pipeline [9]. South Africa was the first African country to adopt the monovalent RVA vaccine (Rotarix®) in September 2009 into its Expanded Program on Immunization (EPI) (WHO, 2009), which culminated in a 77% reduction in RVA disease during the first year that the vaccine was introduced [10,11].

Rotaviruses belong to the *Reoviridae* family. The RV genome is composed of 11 segments of double-stranded RNA (dsRNA) encapsulated in a three-layered protein capsid. Six structural proteins (VP1–VP4, VP6, and VP7) and five or sometimes six non-structural proteins (NSP1–NSP5/NSP6) that encode the RV genome [2]. The outer capsid proteins, VP7 and VP4, which act as neutralizing agents, are universally applied in the binary classification of RV strains into G and P types, respectively [2]. The contemporary classification of RVA strains is based on whole-genome composition underpinned by the nucleotide homology cutoff values that have been determined for the open reading frame (ORF) of each gene segment [12,13]. The numbers of currently described genotypes are 36 G (VP7), 51 P (VP4), 26 I (VP6), 22 R (VP1), 20 C (VP2), 20 M (VP3), 31 A (NSP1), 22 N (NSP2), 22T (NSP3), 27 E (NSP4), and 22 H (NSP5) (http://rega.kuleuven.be/cev/viralmetagenomics/virus-classification).

The globally predominant RVA genotypes are G1P[8], G2P[4], G3P[8], G4P[8], G9P[8], and G12P[8] [14]. However, RVA strains variability by region is well documented [15]. In Africa, RVA genotypes such as G1P[6], G8P[4], G8P[6], G8P[8], and G9P[6] are substantially prevalent but uncommon elsewhere [14–17]. Additionally, G3P[8] and G4P[8] genotypes have been on the decline in Africa and have not been detected in many African countries for almost a decade aside from an impromptu emergence of equine-like G3P[6] and G3P[8] in Botswana and Eswatini [18]. RVAs are classified further into three genogroups: Wa-like, which bears a genotype 1 constellation (I1-R1-C1-M1-A1-N1-T1-E1-H1), DS-1-like, which bears genotype 2 constellation (I2-R2-C2-M2-A2-N2-T2-E2-H2), and a relatively minor AU-1-like characterized by genotype 3 constellation (I3-R3-C3-M3-A3-N3-T3-E3-H3) [19]. Typically, G1P[8], G3P[8], G4P[8], G9P[8], and G12P[8] RVA have a Wa-like genotype constellation, whereas G2P[4], G8P[4], and G8P[6] strains usually have a DS-1-like genotype constellation [19]. G1P[8] is the world's most prevalent genotype accountable for an estimated 50% of RVA infections [20]. The vast antigenic and genetic heterogeneity of G1P[8] strains contributes to the persistent recurrence of VP4 and VP7 protein variants, and the epidemiological fitness of some of these variants might be accountable for their global prevalence [21].

The segmented RNA genome of RVA facilitates reassortment and recombination events, and the error-prone RNA-dependent RNA polymerase promotes high mutation rates [2]. These evolutionary mechanisms lead to the emergence of novel strains and distinct lineages [22]. Intergenogroup reassortment of G1P[8] gene segments has been reported in Africa, Asia, and the Americas [23–27]. These atypical G1P[8] strains were first reported in Okayama Prefecture, Japan during 2012–2013 post-RVA vaccine surveillance of acute gastroenteritis and then in other prefectures, including Aichi, Akita, Kyoto, and Osaka [25–27]. Subsequent incidences were then reported during 2013 post-RVA vaccine surveillance in Phetchabun and Sukhothai provinces in Thailand [28,29] and in 2012–2013

during the pre-RVA vaccine period in Hanoi, Vietnam [30]. Although unpublished, sequence data of G1P[8] DS-1-like sequence strains isolated during pre-vaccine period between August–November 2012 in Palawan, Southwestern region of Philippines have been deposited in the GenBank database. Recently, for the first time in the Americas, G1P[8] DS-1-like strains were reported in 2013 during post-RVA vaccination period from the states of Sao Paulo and Goias in Brazil [23]. In Africa, Jere and colleagues reported the emergence of atypical G1P[8] strains during the post-RVA vaccination period in Blantyre, Malawi [24]. It is not definitively resolved whether these atypical G1P[8] strains are widespread. In addition, there is a paucity of information on whole-genome sequences of G1P[8] strains post-vaccine era with only a few countries performing full-genome characterization of the strain [15,21,31–35]. The African Enteric Viruses Genome Initiative (AEVGI) is conducting whole-genome characterization of country-specific pre- and post-vaccine RVA strains in Africa and has identified, for the first time in South Africa, atypical G1P[8] strains that were circulating before vaccine introduction. This study aimed to determine the genetic relationship and the evolutionary origin of these pre-vaccine atypical G1P[8] RVA strains.

#### **2. Results**

#### *2.1. Nucleotide Sequencing*

Illumina® MiSeq sequencing yielded 14.7 × 10<sup>5</sup> reads (379 bp fragment size) and 11.3 × 10<sup>5</sup> reads (364 bp fragment size) for strains RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/G1P[8], respectively. All the sequences had a phred score of Q ≥ 30 (99.9% base calling accuracy).

#### *2.2. Full-Genome Constellation Analysis*

Whole-gene sequences of the 11 genes of strains, RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/ G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/G1P[8], were determined, and their genotype constellations were revealed as G1-P[8]-I2-R2-C2-M2-A2-N2-T2-E2-H2 (Table 1). The sizes of full-length segments 1 to 11 and their respective open reading frames (ORFs) for the two study strains were determined (Table 1). The ORF sequences for all the 11 genes of these two South African atypical G1P[8] strains were deposited in GenBank under accession numbers MT163245-MT163266.



Color codes indicate genogroup attribution. Green color represents the genotype associated with the Wa-like genogroup, while red color represents the genotype belonging to the DS-1-like genogroup. The nomenclature of the RV strains indicates RV group, species where the strain was isolated, name of the country where the strain was originally isolated, common name, year of isolation, and genotypes for genome segments four and nine as proposed by the Rotavirus Classification Working Group (RCWG) [12]. ORF = open reading frame.

#### *2.3. Sequence and Phylogenetic Analysis*

#### 2.3.1. Phylogenetic Analysis of VP7

Phylogenetically, the diversity of the VP7 G1 genes has been established through seven known lineages (I-VII) [36] (Figure 1). The VP7 genes of the atypical G1P[8] study strains, RVA/ Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/2008/G1P[8], clustered in genetic lineage I, which consisted of a global collection of G1

strains that circulated from 2002 to 2015 (Figure 1). In this lineage I, the two G1 study strains clustered closely together and shared almost absolute gene identities amongst themselves—nt (aa) 99.9% (100%) (Figure 1; Supplementary data 1 (S1)). Analysis of the G1 study strains with locally circulating South African strains retrieved from the GenBank identified the highest sequence identities of 99.8–99.9% (100%) with strain RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P8 and clustered closely with this strain that was isolated the same year, 2008, as the two study strains (Figure 1). However, within the same lineage, the VP7 genes of the two G1 study strains from South Africa clustered distinctly away from the atypical G1P[8] strains reported in Brazil, Japan, Malawi, Philippines, Thailand, and Vietnam [22–26] and displayed overall nt (aa) similarities that ranged from 87.0–98.5% (87.1–98.8%) (Table S1 in Supplementary data 2 (S2)). Specifically, the nt (aa) similarities ranged from 96.9–97.0% (98.5%), 96.9–97.1% (98.5–98.8%), 87.0–98%5 (87.1–98.8), 96.5–96.8% (98.2–98.5%), 96.9–97.0% (98.8%), and 97.0–97.1% (96.9%) to the post-vaccination G1 atypical strains reported in Brazil, Japan, Philippines, Malawi, Thailand, and Vietnam, respectively (S1).

When the two South African G1 study strains were compared to the typical G1 strains selected globally, they displayed the highest nt (aa) similarities of 99.7–99.8% (100%) with a European strain, RVA/Human-wt/BEL/BEL00017/2006/G1P[8] (S1). The nt(aa) similarities comparison to representative strains from Africa (Eastern Africa, Southern Africa, and West Africa), America, Asia, Europe, and Oceania ranged from 93.9–97.6% (94.2–98.2%), 97.1–99.5% (98.2–100%), 97.1–97.9% (96.9–98.2%), 93.1.5–97.6% (94.8%–98.5%), 96.4–99.6% (97.8–99.7%), 99.7–99.8% (100%), and 93.7–98.4% (94.5–98.8%), respectively (Table S2 in S2). In addition, comparison of the VP7 genes of the two study strains to cognate gene sequence of the Rotarix® and RotaTeq® RV vaccine strains displayed nt (aa) identities that ranged from 94.2–94.3% (95.7%) and 91.0–91.1% (93.2%), respectively (S1).

#### Analysis of the VP7 Neutralization Epitopes

The VP7 genes contain three established neutralization epitopes: 7-1a, 7-1b, and 7-2. Twenty-nine amino acids (14 residues in 7-1a, 6 residues in 7-1b, and 9 residues in 7-2) define the three VP7 antigenic epitopes [37]. The VP7 neutralization epitope sites of the two South African study strains were aligned and mapped against cognate neutralization sites of the two RV vaccines, Rotarix®, and RotaTeq®. Four amino acid differences (N94S, S123N, K291R, and M217T) in the VP7 genes of the two South African study strains were identified relative to Rotarix® VP7 neutralization sites, while five amino acid differences (D97E, S123N, K291R, S147N, and M217T) were identified with comparison to RotaTeq® G1 antigenic sites (Figure 2). Antigenically, similar amino acid residues in the VP7 epitopes of the study strains were observed in the corresponding VP7 epitopes of the multiple atypical G1P[8] strains (Figure 2). The VP7 epitopes of the two South African G1 strains were contrasted with those of globally selected lineage I G1 strains. The analysis showed ten amino acid differences (T91N, S94N, D100N, D100E, N123S, R291K, T242A, N147D, L148F, and T217M) (Figure 2).

0.05

● ≥ **Figure 1.** VP7 phylogenetic tree based on the full-length nucleotide sequences. Strains group A rotavirus (RVA)/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/ UFS-NGS-MRC-DPRU1973/2008/G1P[8] are identified by the black filled circular dots (•). Unusual G1P[8] strains from Malawi, Japan, Thailand, Vietnam, Brazil, and Philippine are indicated. Bootstrap values ≥ 70% are shown adjacent to each branch node. Each scale bar indicates the number of nucleotide substitutions per site.


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

#### 2.3.2. Phylogenetic analysis of VP4

The VP4 genes of the two atypical G1P[8] study strains were phylogenetically compared to the four established lineages (l-IV) of the P[8] genotypes [38] (Figure 3). The P[8] genes of the South African strains, RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Humanwt/ZAF/UFS-NGS-MRC-DPRU1973/2008/G1P[8], clustered in lineage III, which consisted of a global collection of P[8] strains that circulated from 2002 to 2014 (Figure 3). Within the P[8]-lineage-III, the two atypical G1P[8] study strains clustered closely together and shared nt (aa) identities of 99.9% (99.7%) amongst themselves (Figure 3; S1). Homology analysis of the P[8] sequences

of the two South African strains with sequences of South African strains retrieved from the GenBank demonstrated the highest nt (aa) sequence identities of 99.6% (99.1–99.4%) with strain RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8] (Figure 3). However, within the same lineage, the VP4 genes of the two atypical strains from South Africa segregated distinctly away from the atypical strains that have been detected in Brazil, Japan, Malawi, Philippines, Thailand, and Vietnam. They exhibited overall nt (aa) similarities that ranged from 95.2–98.0% (95.1–98.6%) (Table S1 in S2). Specifically, the nt (aa) similarities ranged from 97.6–97.8% (98.1–98.6), 98.2–98.5% (98.5%), 95.2–98.0% (95.1–98.2%), 97.7–97.8% (98.1–98.5%), 97.8–98.0% (97.7–98.2%), and 98.1–98.2% (97.9–98.3%) to the post–vaccination atypical strains reported in Brazil, Japan, Malawi, Philippines, Thailand, and Vietnam, respectively (S1). A comparison of the South African P[8] study strains characterized in this study with a global collection of P[8] strains showed their closeness, and the study strains shared the highest nt (aa) similarity of 99.5% (99.1–99.4%) to a Belgian strain, RVA/Human–wt/BEL/BEL00017/2006/G1P[8] (S1). Overall, the nt (aa) similarities in comparison to representative strains from Africa (Eastern Africa, Southern Africa, Western Africa), America, Asia, Europe, and Oceania ranged from 97.4–97.8% (97.8–98.2%), 86.7–99.1% (91.1–99.2%), 86.7–99.1% (91.1–99.2%), 86.6–99.4% (91.0–99.2%), 86.5.–98.9% (91.1–98.7%), 86.8–99.5% (91.2–99.5%), and 86.5–98.5% (91.0%–98.6%), respectively. In addition, the comparison of the atypical VP4 genes to the P[8] genes of the Rotarix® and RotaTeq® vaccine strains displayed nt (aa) identities that ranged from 90.3–90.4% (93.9–94.2%) and 92.3% (95.2%), respectively (S1).

#### Analysis of the VP4 Neutralization Epitopes

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

0.05

● ≥ **Figure 3.** VP4 phylogenetic tree based on the full-length nucleotide sequences. Strains RVA/Humanwt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/ 2008/G1P[8] are identified by the black filled circular dots (•). Unusual G1P[8] strains from Malawi, Japan, Thailand, Vietnam, Brazil, and Philippines are indicated. Bootstrap values ≥ 70% are shown adjacent to each branch node. Each scale bar indicates the number of nucleotide substitutions per site.

**Figure 4.** Alignment of antigenic residues in VP4 between the P[8] component of Rotarix® and RotaTeq® vaccines and wild type P[8] strains. Antigenic residues are divided in four antigenic epitopes in VP8\* (8-1, 8-2, 8-3, and 8-4) and five antigenic epitopes in VP5\* (5-1, 5-2, 5-3, 5-4, and 5-5). Amino acid changes that have been shown to escape neutralization with monoclonal antibodies are indicated with a black dot. Amino acids that differ between Rotarix® and RotaTeq® are indicated in boldface. Green colored residues are residues that are different from Rotarix®, brown colored residues are different from RotaTeq®, and residues colored in sky blue are different from both Rotarix® and RotaTeq®. Dashes (-) indicate no amino acid sequence.

#### 2.3.3. Phylogenetic Analysis of VP1–VP3 and VP6

The evolutionary relationship of the VP1–VP3 and VP6 genes of the two South African study strains with a selection of global RVA strains was performed. The VP1–VP3 and VP6 genes of the two South African study strains clustered closely and displayed nearly absolute gene identities (≥99.9%) amongst each other (Figures S1–S4 in supplementary data 3 (S3)). The VP1–VP3 and VP6 genes of the two South African study strains clustered closely in sublineage composed mainly of locally circulating South African DS-1-like strains and were all found to cluster closely with cognate genes of strain RVA/Human-wt/ZAF/MRC-DPRU2344/2008/G2P[6], which was co-circulating in the population in the same year, 2008, as the study strains (Figures S1–S4 in S3). The VP1–VP3 and VP6 genes of the two study strains were closely related to cognate genes of strain RVA/Humanwt/ZAF/MRC-DPRU2344/2008/G2P[6] with nt (aa) identities ranging from ≥ 99.8–99.9% (≥99.9–100%) for VP1–VP3 and VP6 genes (S1). Phylogenetic relationship of the VP1–VP3 and VP6 genes of the atypical study strains with the cognate genes of the atypical G1P[8] strains reported in Brazil, Japan, Philippines, Thailand, Vietnam, and Malawi exhibited distinct clustering albeit belonging within the same lineage and displayed overall nt (aa) identities that ranged from ≥ 94.4–98.0% (≥ 97.1–100%) (Figure S1C–F in S3; Table S1 in S2). When VP1–VP3 and VP6 genes of the study strains were compared with selected global strains, highest genetic similarities ranging from ≥ 99.4–100% were identified with cognate genes of G3P[6] strains: RVA/Human-RVA/Human-wt/CMR/ES293/2011/G3P[6],

#### RVA/Human-wt/TGO/MRC-DPRU2206/2009/G3G9P[6], RVA/Human-wt/BEL/F01498/2009/G3P[6], and RVA/Human-wt/UGA/MUL-13-166/2013/G3P[6] for VP1–VP3 and VP6, respectively (S1).

#### 2.3.4. Phlylogenetic Analysis of NSP1–NSP5

The NSP1–NSP5 genes of the two South African study strains were highly identical amongst each other with nt (aa) identity value of ≥99.8%, and clustered closely (Figures S5–S9). Close clustering with cognate genes of a locally circulating strain, RVA/Human-wt/ZAF/MRC-DPRU2344/2008/G2P[6], was observed for the all the NSP1–NSP5 genes, and they shared the highest nt (aa) similarities that ranged from ≥ 99.5–100% (99.4–100%) (S1). In contrast, the NSP1–NSP5 genes of the atypical study strains grouped distinctly away from cognate genes of atypical strains reported in Brazil, Japan, Philippines, Thailand, Vietnam, and Malawi and shared an overall nt (aa) similarities that ranged from ≥ 89.0–99.8% (≥94.3–100%) (Figures S5–S9 in S3; Table S1 in S2). Comparison of NSP1–NSP5 genes of the two South African study strains with corresponding selected reference strains collected globally demonstrated nt (aa) identities in the range of 98.6–100% (99.4–100%) with cognate A2, N2, T2, E2, and H2 genes of strains: RVA/Human-wt/BEL/F01498/2009/G3P[6], RVA/Human-wt/KEN/KDH1968/2014/G3P[6], RVA/Human-wt/GHA/GH018-08/2008/G8P[6], RVA/Human-wt/ZMB/MRC-DPRU1673/2009/G2P[4], and RVA/Human-wt/USA/2007769964/2007/G2P[4], respectively (S1).

#### *2.4. Reassortment Analysis*

The concatenated genomes of strains RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/2008/G1P[8] were compared with two South African strains RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8] and RVA/Human-wt/ZAF/MRC-DPRU2344/G2P[6] (Figure 5). The two atypical South African study strains shared a highly conserved backbone with all genes exhibiting > 99.8% nucleotide similarity. The VP7 and the VP4 genes of the two atypical strains shared the highest genetic similarities to RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8]. However, the internal backbone genes were extremely diverse. The internal backbone genes of the atypical strains exhibited highest genetic similarity to RVA/Human-wt/ZAF/ MRC-DPRU2344/G2P[6]. The results of this analysis suggest that the atypical G1P[8] strains were likely derived via reassortment events between contemporary, endemic South African strains.

**Figure 5.** Nucleotide sequence similarities of the concatenated genome of RVA//Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] were compared with South African strains RVA//Human-wt/ ZAF/UFS-NGS-MRC-DPRU1973/2008/G1P[8], RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8], and RVA/Human-wt/ZAF/MRC-DPRU2344/G2P[6]. The left axis displays the strains, and rotavirus genome segment is included in the top scale. The bottom scale shows distance in kb.

#### **3. Discussion**

This study described the first pre-vaccine era atypical reassortant G1P[8] strains, whose outer gene segments (VP7 and VP4) expressed a Wa-like genotype, whereas the backbone genes expressed a DS-1-like genotype constellation. Analysis of the whole-genome constellation showed that genetic reassortment mechanism generated the DS-1-like G1P[8] strains. Rotavirus reassortment events are mainly facilitated by the segmentation inherent in the RV genome [2], which can generate rare or novel RV strains and hence contribute to the vast RVA diversity [22]. Wa-like and DS-1-like intergenogroup reassortment events involving G1P[8] and DS-1-like genotype constellation have been described recently in six countries: Brazil, Japan, Philippines, Thailand, Vietnam, and Malawi [23–30]. According to literature, viable atypical reassortant strains can occur under natural conditions involving Wa-like G1P[8] or G3P[8] outer capsid genes expressing DS-1-like genetic background [41–44]. This study identified two DS-1-like G1P[8] strains in the course of the ongoing AEVGI whole-genome characterization of South African RVA strains. While reported in low frequencies and limited settings in Brazil (1.6% during 2013–2017 seasons) [22], Thailand (0.4% during 2012–2014 seasons) [27], and Vietnam (14% during 2012/2013 season) [28,29], 31–62% of these DS-1-like G1P[8] strains accounted for RVA positive strains circulating across selected regions in Japan [26], and 40% of randomly sampled post-vaccine samples were reported in Malawi [24]. Such atypical reassortant strains have the potential to predominate in circulation. G1P[4] strains suggested to have emerged from intergenogroup reassortment events accounted for 41% of RVA strains circulating in the peak months of the 2001 RVA season in Detroit, USA [45], whereas a surge in G3P[4] strains also presumed to have emerged from intergenogroup reassortment events were detected in Brazil at 36% [46] and in Ghana at 64% [47]. The two South African study strains were identified during the pre-RVA vaccination period in South Africa in contrast to the previously reported atypical strains found during the post-RVA vaccination period. This implies that the reassortment events that led to the emergence of the South African atypical G1P[8] strains may not necessarily be driven by vaccine-induced selective pressure but by natural evolutionary processes of RVA genome.

In order to identify the ancestral origin of these G1P[8] strains, assessment of whole-gene sequences and phylogenetic analysis showed that the outer capsid genes, VP7 and VP4, of strains RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/2008/G1P[8] were 99.6–99.9% (99.1–100%) identical with South African strain RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8] and clustered together in the same clade within the same. For the internal genes, the highest nucleotide identities were identified with cognate genes of another South African strain, RVA/Human-wt/ZAF/MRC-DPRU2344/2008/G2P[6], detected in the 2008 RVA season. Put together, it is probable that a locally circulating G2P[6] strain such as RVA/Human-wt/ZAF/MRC-DPRU2344/2008/G2P[6] with a DS-1-like backbone derived the VP7 and the VP4 genes from a locally co-circulating G1P[8] strain such as RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8], generating the double-gene reassortants. Consequently, the results obtained in this study indicate that the two atypical South African G1P[8] strains were generated locally through genetic reassortment events. The generation of these reassortant double-gene strains tend to be independent of the events in which Brazilian, Japanese, Thai, Vietnamese, and Malawian DS-1-like G1P[8] strains were generated. The nine internal genes of the South African DS-1-like G1P[8] strains always clustered together with cognate genes of a locally circulating G2P[6] strain distinctly away from the cluster comprising Japanese, Thai, Philippines, and Malawi DS-1-like G1P[8] strains. Therefore, the South African DS-1-like G1P[8] strains emerged clonally from independent events, a phenomenon observed for the Malawian [24] and the Vietnamese [30] DS-1-like G1P[8] strains. In contrast, Brazilian, Japanese, and Thai G1P[8] DS-1-like strains were established to have been derived from a common ancestor [23,29].

Vaccine escape mutants can result due to mutations occurring in well-known VP7 neutralization epitope regions [48]. Host–antigen binding interactions involving human G1 strains are significantly impacted by mutations occurring at positions 94, 97, 147, and 291 [48]. The identified N94S substitution involving the substitution of asparagine (N) with a serine(S), which are both polar non-charged amino acid residues [49], may not significantly alter the overall morphology of the protein surface. However, since asparagine is usually N-glycosylated, there is a likely loss of glycosylation site, which could have a wide-ranging impact on the immunogenicity of the 7-1a epitope [50]. The D97E amino acid substitution involving polar negatively charged residues, aspartate (D) and glutamate (E), is likely to be a silent nucleotide change [49]. Similarly, a K291R substitution involving lysine (K), an amphipathic

polar amino acid, and arginine (R), a positively charged amino acid, is unlikely to have a far-reaching structural effect on the VP7 protein surfaces. However, M217T substitution was identified resulting in substitution of methionine, a non-polar residue to threonine, a polar amino acid residue, which could likely result in significant changes in biochemical properties of VP7 [49]. This M217T substitution was also present in earlier strains as well as some post-vaccine strains that were included for analysis, and the role it plays in driving epidemiological fitness of G1 strains is not fully resolved. In the host cell, trypsin-like proteases cleave the VP4 spike protein into two structural domains (VP8\* and VP5\*) [2]. Four surface-exposed antigenic epitopes (8-1 to 8-4) have been described in the VP8\* region, while five antigenic epitopes (5-1 to 5-5) in the VP5\* region have been documented [40]. The amino acid changes E150D, N195G, S125N, and N135D that were observed relative to the vaccine strains were conservative. However, a S131R substitution that resulted in a change in polarity might play a role in escape of host immunity [49]. Another amino acid substitution R131S resulting in a change in charge from positively charged amino acid to non-charged amino acid that was identified when comparison was made against globally selected lineage-III VP4 strains might impact vaccine escape effect [49].

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

#### *4.1. Ethics Approval*

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

#### *4.2. Sample Collection*

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

#### *4.3. Extraction and Purification of Double-Stranded RNA*

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

integrity and the enrichment of dsRNA, which was visualized on a G:Box Syngene UV transilluminator (Syngene, Cambridge, UK).

#### *4.4. Synthesis and Purification of Complementary DNA (cDNA)*

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

#### *4.5. DNA Library Preparation and Whole-Genome Sequencing*

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

#### *4.6. Genome Assembly*

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

#### *4.7. Determination of Rotavirus Whole-Genotype Constellations*

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

#### *4.8. Phylogenetic Analyses*

Complete sequences for each gene segment were aligned and sequence comparisons performed as described previously [53–55]. Multiple sequence alignments were implemented utilizing the MUSCLE package in Molecular Evolutionary Genetics Analysis (MEGA) 6 software ([56]; http: //www.megasoftware.net/). Upon alignment, the DNA Model Test program in MEGA 6 was used to determine the evolutionary model that best fits each gene sequence datasets. The models identified as best fitting with the sequence data for the indicated genes using the Corrected Akaike Information Criterion (AICc) were as follows: GTR+G+I (VP7, VP4, VP6, VP1, VP2, VP3, NSP1, NSP2, NSP3) and

HKY+G+I (NSP4 and NSP5). These models were utilized in maximum-likelihood trees' construction using MEGA 6 with 1000 bootstrap replicates to estimate branch support. Genetic distance matrices were prepared using the *p*-distance algorithm of MEGA 6 software [56]. In addition to the two whole-genome sequences of the strains in this study, other cognate sequences were acquired from GenBank ([57]; http://www.ncbi.nlm.gov/genbank). Further phylogenetic analysis by geographical regions Africa (Eastern Africa, Southern Africa, and West Africa), Asia, Americas, Europe, and Oceania was also performed. mVISTA software was used to visualize the comparative sequence similarities of concatenated whole-genome of genetically related strains [58].

#### **5. Conclusions**

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

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/5/391/s1, Supplementary data 1 (S1): Identity matrices analysis for VP1, VP2, VP3, VP4, VP6, NSP1, NSP2, NSP3, NSP5 and NSP5 nucleotide and deduced amino acid identities among strains calculated by distance matrices using P-distance algorithm in MEGA 6. Supplementary data 2 (S2): Homology analysis summary containing Table S1 and Table S2. Supplementary data 3 (S3): Additional phylograms containing Figure S1:VP1, S2:VP2, S3:VP3, S4:VP6, S5:NSP1, S6:NSP2, S7:NSP3, S8:NSP4 and S9:NSP5.

**Author Contributions:** M.M.N., K.C.J., F.E.D. and V.N.N. conceptualized the main project. P.N.M., M.T.M., M.M.N. and S.P.R. performed the laboratory experiments. M.J.M. and M.L.S. facilitated the sample resources. Formal analysis was done by P.N.M. and M.M.N. Data curation was performed by P.N.M., M.T.M., M.M.N. and S.P.R. Writing of the original draft preparation was performed by P.N.M. Review of the drafts was performed by all co-authors. Supervision, funding acquisition and project administration was performed by M.M.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was principally funded by a grant awarded to M.M.N. by Bill and Melinda Gates Foundation (BMGF-OPP1180423\_2017). Other funding grants awarded to M.M.N. that funded this research include: South African Medical Research Council (SAMRC) through the Self-Initiated Research grant (SIR), Poliomyelitis Research Foundation (PRF-19/16) and National Research Foundation (NRF-120814).

**Acknowledgments:** Assistance in retrieval of the archived stool samples by Khutso Mothapo, Kebareng Rakau and Nonkululeko Magagula at the WHO-RRL in Pretoria, is hereby acknowledged. Assistance in laboratory work by Lesedi Mosime, Gilmore Pambuka and Teboho Mooko is duly acknowledged. Some aspects of data analysis guidance by Mathew Esona (CDC, Atlanta) and Felicity Burt is greatly acknowledged. Technical ICT support by Stephanus Riekert is also acknowledged.

**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*

### **Molecular Characterisation of a Rare Reassortant Porcine-Like G5P[6] Rotavirus Strain Detected in an Unvaccinated Child in Kasama, Zambia**

**Wairimu M. Maringa <sup>1</sup> , Peter N. Mwangi <sup>1</sup> , Julia Simwaka <sup>2</sup> , Evans M. Mpabalwani <sup>3</sup> , Jason M. Mwenda <sup>4</sup> , Ina Peenze <sup>5</sup> , Mathew D. Esona <sup>5</sup> , M. Je**ff**rey Mphahlele 5,6 , Mapaseka L. Seheri <sup>5</sup> and Martin M. Nyaga 1,\***


Received: 24 July 2020; Accepted: 14 August 2020; Published: 17 August 2020

**Abstract:** A human-porcine reassortant strain, RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/ G5P[6], was identified in a sample collected in 2014 from an unvaccinated 12 month old male hospitalised for gastroenteritis in Zambia. We sequenced and characterised the complete genome of this strain which presented the constellation: G5-P[6]-I1-R1-C1-M1-A8-N1-T1-E1-H1. The genotype A8 is often observed in porcine strains. Phylogenetic analyses showed that VP6, VP7, NSP2, NSP4, and NSP5 genes were closely related to cognate gene sequences of porcine strains (e.g., RVA/Pig-wt/CHN/DZ-2/2013/G5P[X] for VP7) from the NCBI database, while VP1, VP3, VP4, and NSP3 were closely related to porcine-like human strains (e.g., RVA/Human-wt/CHN/E931/2008/G4P[6] for VP1, and VP3). On the other hand, the origin of the VP2 was not clear from our analyses, as it was not only close to both porcine (e.g., RVA/Pig-tc/CHN/SWU-1C/2018/G9P[13]) and porcine-like human strains (e.g., RVA/Human-wt/LKA/R1207/2009/G4P[6]) but also to three human strains (e.g., RVA/Human-wt/USA/1476/1974/G1P[8]). The VP7 gene was located in lineage II that comprised only porcine strains, which suggests the occurrence of independent porcine-to-human reassortment events. The study strain may have collectively been derived through interspecies transmission, or through reassortment event(s) involving strains of porcine and porcine-like human origin. The results of this study underline the importance of whole-genome characterisation of rotavirus strains and provide insights into interspecies transmissions from porcine to humans.

**Keywords:** whole-genome; genotype constellation; interspecies transmission; reassortment; porcine; porcine-like human

#### **1. Introduction**

Group A rotaviruses (RVA), of the family *Reoviridae*, are the number one viral pathogens causing severe diarrhoea in children below five years of age [1]. In 2016, an estimated 128,000 deaths in children below five years were due to RVA infections, 90% of which occurred in developing countries [2,3]. Similarly, RVA are the primary cause of acute gastroenteritis in new-born piglets [4].

Rotaviruses have a distinctive morphology which comprises a nonenveloped, three-layered icosahedral protein shell. The rotavirus genome within the protein shell comprises 11 segments of double-stranded (dsRNA) that encode six structural viral proteins (VP1 to VP4, VP6, and VP7) and five or six nonstructural proteins (NSP1 to NSP5/6) [1]. A binary classification system is used to distinguish RVA based on the antigenic properties of the outer shell proteins, VP7 and VP4, that determine the G-genotype and P-genotype, respectively [1]. Furthermore, RVA can be separated into two main genogroups and one minor genogroup according to a whole-genome classification system, whereby a specific genotype is assigned to the 11 gene segments. These genogroups represent the genotype constellations that are present in most human strains globally [5,6]. Genogroup 1 (Wa-like) bears the constellation I1-R1-C1-M1-A1-N1-T1-E1-H1 and is often associated with the G genotypes G1, G3, G4, G9, and G12 and P genotype P[8]. Genogroup 2 (DS-1-like) includes G2P[4] strains and bears the constellation I2-R2-C2-M2-A2-N2-T2-E2-H2. Lastly, the minor genogroup 3 (AU-1-like) bears the I3-R3-C3-M3-A3-N3-T3-E3-H3 constellation and includes G3P[9] strains [7]. As of 5th May 2020, the Rotavirus Classification Working Group had identified at least 36 G, 51 P, 26 I, 22 R, 20 C, 20 M, 31 A, 22 N, 22 T, 27 E, and 22 H genotypes [8]. The whole-genome classification system has made it possible to analyse and understand the origin of various strains, interspecies transmission, and animal–human reassortment events [9]. Human Wa-like strains and porcine rotavirus strains share a common origin, whereas DS-1-like and AU-1-like strains have a common origin with bovine and feline strains, respectively [5].

In humans, G1-G4, G9, and G12 along with P[4], P[6], and P[8] are the most frequently detected, globally [10–13]. On the contrary, in porcine, predominant genotypes are G3-G5, G9, and G11 along with P[6], P[7], and P[13] [4,14]. Porcine rotaviruses bear the constellation I5-R1-C1-M1- A8-N1-T1/T7-E1-H1 [5,15–20]. While human Wa-like RVA differ from porcine rotaviruses in some gene segments (VP4, VP6, VP7, and NSP1), they both appear to have genotype 1 in the VP1, VP2, VP3, NSP2, NSP3, NSP4, and NSP5 gene segments. Hence, the suggestion that human Wa-like and porcine RVAs have arisen from a common ancestor [5].

The findings that show animals can serve as potential reservoirs for genetically diverse rotavirus strains that can be passed on to humans have elicited a large amount of interest and topics for further research [21]. Several novel and rare animal-like or animal–human reassortant rotavirus strains have been identified globally [22–28]. The detection of animal strains in humans is presumed to be as a result of zoonotic transmission, along with reassortment, which contributes to the diversity of circulating RVA [4,29,30]. Inter- and intragenogroup reassortment may occur when multiple RVA simultaneously infect a host. This is attributed to the segmented nature of the rotavirus genome [1,31]. It is, therefore, necessary to continuously carry out the monitoring of animal RVA and the role they play in contributing to the diversity of circulating RVA in humans.

The G5, one of the most common porcine genotypes, has sporadically been identified in human populations in Brazil (G5P[X]), Cameroon (G5P[7] and G5P[8]), Argentina (G5P[8]), and the United Kingdom(G5P[X]) [32–36]. The P[6] is presumed to be of porcine origin. They have also been identified in human populations [37–40]. The first human G5P[6] strain, LL36755, was detected in a child who had acute gastroenteritis in China in 2007 [41]. Other G5P[6] strains were detected in Vietnam, Taiwan, Bulgaria, Japan, and Thailand [37,42–45]. To date, the whole-genome of only two human G5P[6] strains—Bulgarian BG620 (nt sequences unavailable in the DDBJ, EMBL, and GenBank data libraries as of 13 August 2020) and Japanese Ryukyu-1120 (full open reading frame, available in GenBank)—have been analysed [45,46].

Diarrhoea is a burden for the Zambian healthcare system, with about 33% of the extreme cases being attributable to RVA [47–49]. In an attempt to generate disease burden attributable to rotavirus diarrhoea in children, the Zambian Ministry of Health, with support from WHO, launched rotavirus surveillance at the University Teaching Hospital (UTH) in 2006 [50,51]. Surveillance data generated provided evidence of the burden of rotavirus diarrhoea that supported the introduction of the rotavirus vaccine, Rotarix®, as a pilot project in Lusaka, Zambia in 2012, and was later rolled out nationwide in November 2013 [50]. According to the estimates reported by the World Health Organization (WHO) and the United Nations International Children's Emergency Fund (WHO/UNICEF), rotavirus vaccine coverage in Zambia has been consistently high for the last six years, increasing from 73% in 2014 to 90% in 2019 [52]. Over this period, a sustained and significant reduction in rotavirus-associated hospitalisations and mortality was observed in children under 5 years [51].

The African Rotavirus Surveillance Network, coordinated by the World Health Organization Regional Office for Africa (WHO/AFRO), is actively monitoring the diversity and distribution of RVA genotypes in children hospitalised with acute diarrhoea [53]. Initially, the network was established with four countries in 2006, and expanded to 29 countries by the end of 2016 [54,55]. The Diarrhoeal Pathogens Research Unit at Sefako Makgatho University in Pretoria (South Africa) and the Noguchi Memorial Institute for Medical Research in Accra (Ghana) are the two WHO Rotavirus Regional Reference Laboratories (RRLs) for the network that conducts monitoring of rotavirus epidemiology in Africa [55]. The WHO/AFRO is currently supporting the University of the Free State-Next Generation Sequencing (UFS-NGS) unit to undertake rotavirus surveillance of rotavirus strains that circulated in Zambia between 2013 and 2016 at the whole-genome level. A G5P[6] strain, UFS-NGS-MRC-DPRU4723, was identified among these strains and was analysed so as to elucidate its origin and evolution. The sample was collected in 2014 from an unvaccinated 12 month old male hospitalised for gastroenteritis at Arthur Davison Children's Hospital in Ndola, Zambia.

#### **2. Results**

#### *2.1. Nucleotide Sequencing and Identity of the Strain*

Illumina® MiSeq sequencing exhibited a phred score of Q30 and collectively yielded 98.8Mbs of data for this specific sample. The whole genome of RVA/Human-wt/ ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] was 18272 bps in size. The length and ORF of the 11 gene segments as determined by nucleotide sequencing are shown in Table 1. A BLASTn search was performed, and it appeared to exhibit maximum sequence identities of 95.7%–98.0% with porcine and human porcine-like strains (Table 1). Based on the whole genome classification system, RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] exhibited a G5-P[6]-I1-R1-C1-M1-A8-N1-T1-E1-H1 genotype constellation (Table 2). The genetic constellation of the study strain was compared to those of other G5 and non-G5 strains retrieved from the GenBank (Table 2).


**Table 1.** The segment and ORF lengths of strain UFS-NGS-MRC-DPRU4723 and the highest sequence identities obtained using the Basic Local Alignment Search Tool (BLAST).

#### *2.2. Sequence and Phylogenetic Analysis*

To investigate the potential origin of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6], phylogenetic trees were constructed for each of the 11 gene segments along with cognate gene sequences of RVA strains obtained from the GenBank.

#### 2.2.1. Sequence and Phylogenetic Analysis of the VP7 Gene

Phylogenetically, there are three known VP7 G5 lineages (I-III) [63]. The VP7 genes of RVA/Human-wt/ ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered into lineage II, which consisted only of porcine G5 strains from mainly Asia and the Americas (Figure 1). The VP7 gene showed the highest nucleotide (nt) and amino acid (aa) identities with the Chinese porcine strains RVA/Pig-wt/CHN/DZ-2/2013/ G5P[X] nt (aa), 98.6% (99.0%), and RVA/Pig-wt/CHN/JN-2/2014/G5P[X] 98.5% (99.0%) and was distantly related to the strains within lineage III with lower sequence identities (nt, 83.4%–86.5%; aa, 90.4%–94.5%) (Figure 1; Supplementary data 1). Overall, strains within lineage II exhibited sequence identities that were in the range nt, 89.6%–98.6%; aa, 92.4%–99.0% (Supplementary data 1).

The comparison of the amino acid sequence of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/ 2014/G5P[6] to reference G5 strains e.g., RVA/Pig-wt/THA/CMP-001-12/2012/G5P[13] (lineage I), RVA/Pig-wt/BRA/ROTA24/2013/G5P[6] (lineage II) and RVA/Human-wt/JPN/Ryukyu-1120/2011/G5P[6] (lineage III) within each of the three lineages revealed a high identity (range 90.0%–94.9% (Supplementary data 1; Supplementary data 2a). Numerous substitutions were identified in the nine VP7 variable regions, VR-1 to VR-9 [64]: VR-1 (I9V and I19V), VR-2 (V27T and V29T), VR-3 M/F39L, I40V, V41I, L/I43V, I/L/V47F, R49K, and A50T), VR-4 (K/A65T, V/M68A, M/A72T, and M/Q75T), VR-5/antigenic site A (N/S/D/T96A), VR-6 (I129V and D130E), VR-7/antigenic site B (N145D and A/V/E146G), VR-8/antigenic site C (L/S208T, A210T, T/V212I, S/A213I, I/M217T, V218I, and S220N), and VR-9/antigenic site F (A/M241T and S242N).


**Table 2.** Genotype natures of the 11 gene segments of Zambian strain UFS-NGS-MRC-DPRU4723 compared with those of selected human and porcine strains.





Blue shading indicates the gene segments with genotypes identical to those of UFS-NGS-MRC-DPRU4723. Bold font indicates genotypes associated with porcine strains. "−" indicates that no sequence data were available in GenBank/EMBL/DDBJ data banks. \* Genotype assignment based on reports by [37] (strain 03-98sP50) and (strain BG260) [46]. To date, the nucleotide accession numbers for the 11 gene segments of strains 03-98sP50 and BG260 are not available in the GenBank, EMBL, or DDBJ data banks.

0.05

‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ▪ ≥ **Figure 1.** Phylogenetic tree constructed from the nucleotide sequences of the VP7 genes of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] and representative strains. The position of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] is shown by the black square (■). Reference strains obtained from GenBank are represented by accession number, strain name, country, and year of isolation. The three closest strains, as identified by BLASTn, are also included. Bootstrap values ≥70% are shown adjacent to each branch node. Scale bar: 0.05 substitutions per nucleotide.

#### 2.2.2. Sequence and Phylogenetic Analysis of the VP4 Gene

The VP4 gene of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] was phylogenetically compared to the already established five lineages (I-V) of genotype P[6] [65] (Figure 2). The P[6] gene of the study strain clustered into lineage V, which consisted of porcine and putative human porcine-like strains detected in parts of Europe and one African strain. A similarity analysis of the P[6] gene of the study strain with strains obtained from GenBank showed that the Zambian G5P[6] exhibited the highest sequence identity of 98.1% (98.3%) with a porcine-like human strain RVA/Human-wt/COD/KisB332/2008/G4P[6] from the Democratic Republic of Congo (Supplementary data 1). All the African strains clustered into a separate lineage, lineage I, with sequence identities of 85.7%–86.8% (92.5%–93.9%) (Supplementary data 1).

The deduced amino acid sequences of the VP4 gene of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] along with the reference P[6] strain from each of the five lineages was compared (Supplementary data 2b). The reference strains shared high amino acid identities ranging from 91.0% to 98.3% (Supplementary data 1). Several amino acid changes were identified throughout the VP4 protein, and most of the substitutions were concentrated in the hypervariable region (amino acid 71-208) which houses the VR-3 (92–192) and includes a neutralization site at amino acid 135 [66,67]. Several amino acid substitutions were observed among the P[6] lineage I strains [65] at the VR-3 (L105I, V108I and T134S) and VR-8 (D602N) variable regions. Other amino acid substitutions were identified among the P[6] lineages at VR-1 (S30N), VR-2 (I61V), VR-3 (V112I, N114S, V130I, H182N and T189S), VR-4 (I280V), and VR-9 (E698K). The potential trypsin cleavage sites at residues 241 and 247 [68] were highly conserved in all the strains with three substitutions at positions 242 (I to V), 243 (A to T), and 244 (H to Y).

#### 2.2.3. Phylogenetic Analysis of the VP6 Gene

The VP6 gene of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered closely with divergent African porcine strains from Uganda (RVA/Pig-wt/UGA/BUW-14-A003/2014/G3P[13], RVA/Pig-wt/UGA/KYE-14-A048/2014/G3P[13], and RVA/Pig-wt/UGA/KYE-14-A047/2014/G3P[13]) and a human porcine-like strain from the Democratic Republic of Congo (RVA/Human-wt/COD/ KisB332/2008/G4P[6]) which displayed nt(aa) sequence identities ranging from 98.6% to 98.9% (98.9%–99.7%) (Figure 3, Supplementary data 1). Porcine-like Asian strains such as RVA/Human-wt/CHN/ GX54/2010/G4P[6] and RVA/Human-wt/CHN/E931/2008/G4P[6] clustered separately, displaying identities of 88.7%–90.2% (97.5%–98.7%) (Supplementary data 1).

‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ▪ ≥ **Figure 2.** Phylogenetic tree constructed from the nucleotide sequences of the VP4 genes of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] and representative strains. The position of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] is shown by the black square (■). Reference strains obtained from GenBank are represented by accession number, strain name, country, and year of isolation. The three closest strains, as identified by BLASTn, are also included. Bootstrap values ≥70% are shown adjacent to each branch node. Scale bar: 0.05 substitutions per nucleotide.

0.05

0.05

‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ▪ ≥ **Figure 3.** Phylogenetic tree constructed from the nucleotide sequences of the VP6 genes of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] and representative strains. The position of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] is shown by the black square (■). Reference strains obtained from GenBank are represented by accession number, strain name, country, and year of isolation. The three closest strains, as identified by BLASTn, are also included. Bootstrap values ≥70% are shown adjacent to each branch node. Scale bar: 0.05 substitutions per nucleotide.

#### 2.2.4. Phylogenetic Analysis of VP1 Gene

The VP1 gene of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered only with porcine and porcine-like human strains from Asia (China and Vietnam) (Supplementary data 3a). The VP1 gene exhibited a maximum nt (aa) sequence identity of 96.8% (98.9%) with the Chinese human porcine-like reassortant strains RVA/Human-wt/CHN/GX82/2010/G4P[6], RVA/Human-wt/CHN/GX78/2010/G4P[6], RVA/Human-wt/CHN/GX77/2010/G4P[6], and RVA/Humanwt/CHN/GX54/2010/G4P[6] (Supplementary data 1). Overall, the Asian strains within the cluster showed sequence identities of 94.1%–96.8% (97.9%–98.9%). Human non-porcine African strains clustered separately, with lower identities of 88.2%–88.8% (96.3%–97.3%) (Supplementary data 1).

#### 2.2.5. Phylogenetic Analysis of VP2 Gene

The VP2 gene of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] fell into a distinct cluster predominantly composed of porcine and porcine-like human strains from Asia (China, India, Vietnam, South Korea, and Sri Lanka) (Supplementary data 3b). The VP2 gene of the study strain showed a maximum nt (aa) sequence identity of 96.6% (90.9%) with a Sri Lankan porcine-like human strain RVA/Human-wt/LKA/R1207/2009/G4P[6] (Supplementary data 1).

#### 2.2.6. Phylogenetic Analysis of VP3 Gene

The VP3 gene of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered in a lineage composed mainly of Asian (Asia and Thailand) porcine and porcine-like human strains (Supplementary data 3c), and exhibited the highest nt (aa) sequence identity with the Chinese porcine-like human strains—RVA/Human-wt/CHN/R946/2006/G3P[6], 95.8% (97.8%) and RVA/Human-wt/CHN/E931/2008/G4P[6], 95.7% (98.0%) (Supplementary data 1). The overall similarities of the Asian strains within the lineage ranged from 84.8% to 95.8% (92.7%–97.8%) (Supplementary data 1). Non-porcine African strains clustered separately and showed lower sequence identities of 84.1%–84.5% (92.1%–92.7%) (Supplementary data 1).

#### 2.2.7. Phylogenetic Analysis of NSP1 Gene

The NSP1 gene of strain RVA/Human-wt/ZMB/UFS-NSG-MRC-DPRU4723/2014/G5P[6] was assigned to a porcine genotype A8 and clustered among Asian (Vietnam, China, and Bangladesh) porcine and porcine-like human strains and an African (Ghana) porcine strain (Supplementary data 3d). The NSP1 gene of the study strain was closest to strain RVA/Human-tc/VNM/NT0042/2007/G4P[6] displaying a nt(aa) sequence identity of 98.2% (97.9%) (Supplementary data 1). The porcine and porcine-like human strains from Europe and the Americas clustered separately showing sequence identities of 84.2%–85.9% (85.4%–88.2%) and 84.1%–85.9% (83.7%–88.3%), respectively (Supplementary data 1).

#### 2.2.8. Phylogenetic Analysis of NSP2 Gene

The NSP2 gene of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered with Asian and European porcine and porcine-like human strains (Supplementary data 3e). The Nt(aa) similarity analysis showed that the NSP2 gene of the study strain was most similar to the Chinese porcine strains RVA/Pig-wt/CHN/YN/2012/GXP[X] and RVA/Pig-tc/CHN/SCMY-A3/2017/G9P[23]—96.8% (97.8%) (Supplementary data 1). Two African porcine strains, RVA/Pig-wt/ZAF/MRC-DPRU1487/2007/G3G5P[23] and RVA/Pig-wt/ZAF/MRC-DPRU1557/2008/G4G5P[23], were seen to cluster within the same lineage with sequence identities of 93.6%–93.7% (97.5%–97.8%) (Supplementary data 1).

#### 2.2.9. Phylogenetic Analysis of NSP3 Gene

The NSP3 gene of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered closely with porcine and porcine-like human strains mainly from Asia (Thailand and Vietnam) and exhibited a maximum nt(aa) sequence identities of 96.5%–97.0% (98.4%–98.7%) with the strains RVA/Human-wt/VNM/30378/2009/G26P[19], RVA/Pig-wt/VNM/12070-4/2012/GXP[X], RVA/Human-wt/VNM/NT0205/2007/G4P[6], and RVA/Human-wt/VNM/NT0621/2008/G4P[6] (Supplementary data 1; Supplementary data 3f). ‐ ‐ ‐ ‐ ‐

‐ ‐ ‐ ‐

#### 2.2.10. Phylogenetic Analysis of NSP4 Gene

The NSP4 gene of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered with porcine and porcine-like human strains identified in Asia (China and Vietnam) and a porcine-like human strain from the Americas (Brazil) (Supplementary data 3g). In this cluster, the closest strains to UFS-NGS-MRC-DPRU4723 were the wild pig strains (RVA/WildBoar-wt/CZE/P828/2015/G9P[23] and RVA/WildBoar-wt/CZE/P830/2015/G9P[23]) from the Czech Republic, with nt(aa) sequence identities of 97.5% (98.3%) (Supplementary data 1). The Asian strains within the cluster showed nt(aa) similarities of 96.2%–97.3% (97.7%–98.9%). Porcine and porcine-like human strains from the Americas clustered separately and exhibited identities of 87.2%–96.4% (94.3%–98.9%) (Supplementary data 1). ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ 

#### 2.2.11. Phylogenetic Analysis of the NSP5 Gene

The NSP5 gene of strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] clustered with porcine strains from Asia and showed the highest nt(aa) sequence identity of 98.6% (100%) with the porcine strains RVA/Pig-wt/CHN/TM-a/2009/G3P[8] and RVA/Pig-tc/CHN/TM-a-P20/2018/G9P[23] identified in China (Supplementary data 1; Supplementary data 3h). Overall, the porcine and porcine-like human strains from Asia and the Americas displayed nt(aa) identities of in the range 94.8%–98.6% (98.0%–100%) and 93.9%–96.1% (95.9%–99.0%), respectively (Supplementary data 1). ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐

#### *2.3. Reassortment Analysis*

The concatenated whole genome alignment of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/ 2014/G5P[6], together with the Japanese G5P[6] strain and selected Chinese porcine-like human P[6] strains, was visualised (Figure 4). The whole genome of the Zambian G5P[6] strain demonstrated a relatively high degree of conservation with the Japanese G5P[6] strain and the two Chinese G4P[6] strains. With the exception of VP7 and VP4, the genome of the Chinese strain E931 exhibited the overall highest genomic conservation to the study strain. With the exception of VP7, VP3, and NSP1 genes, the Chinese strain GX54 shared a highly conserved genome with the study strain. The Japanese strain Ryukyu-1120 demonstrated a highly similar genome to the study strain for seven of the 11 genes, the exceptions being VP1, VP3, VP6, and VP7. The results of this analysis confirmed the genetic similarity between RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] and Asian (Chinese) porcine-like human strains, hence suggesting that the Zambian G5P[6] strain may have been derived via reassortment events. ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐

‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ ‐ **Figure 4.** mVISTA whole genome nucleotide alignment comparing the Zambian G5P[6] strain (RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014G5P[6]) with the G5P[6] strain from Japan (Ryukyu-1120), whose whole genome sequence had been determined, and with selected porcine-like human P[6] strains from China (GX54 and E931). Strain names are shown on the left, and the proteins VP1-VP4, VP6-VP7, and NSP1-NSP5 are indicated on the top. The bottom scale indicates distance in kb. Percentile values on the right indicate sequence-based similarity between the study strain and the respective reference strains. Shading indicates the level of conservation.

#### **3. Discussion**

The detection of genotype G5 in humans, which is typical for pigs, is possibly due to interspecies transmission [35,45]. In Zambia, as with many countries in Africa, humans and farm animals live in proximity. The interaction between humans and animals could be the primary cause for zoonotic transmission, which could result in genetic reassortments and perhaps other mechanisms of genetic diversity, ultimately leading to the introduction and spread of animal genotypes into human populations [69].

In this study, an analysis was conducted on a sample collected from a child admitted to a paediatric ward presenting with clinical symptoms (vomiting, diarrhoea, and fever) that are usually present during typical rotavirus infection. This raises the question whether such animal-derived strains are capable of mutating and effectively spreading within/across human populations as in the case of established typical Wa-like and DS-1-like genotype constellations, with the same magnitude of rotavirus disease severity. Furthermore, taking into consideration that the G5 and P[6] genotypes are not included in the currently available vaccines, the probability for such strains to have the potential to spread more swiftly from human to human may have implications for the effectiveness of current rotavirus vaccine candidates that are in use in African countries.

This study identified the complete genome of a reassortant porcine-like human strain, G5P[6], that showed the genotype constellation G5-P[6]-I1-R1-C1-M1-A8-N1-T1-E1-H1, which is commonly found in porcine and porcine-like human rotavirus strains [19]. RVA/Human-wt/ ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] was found to share the same constellation (I1-R1-C1- M1-A8-N1-T1-E1-H1) with the archival porcine strain, Gottfried, and porcine-like human strains—BG260, E931, and GX54 [5,46,56,58]. In addition, porcine strains 12R002, 12R005, and 12R006, as well as porcine-like human strains Ryukyu-1120, mani-97, 30378, rj24598, and BE2001 shared the same constellation with strain RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] with the exception of VP6 (I5 instead of I1) and NSP3 (T7 instead of T1 gene segments) [20,25,26,45,70].

A phylogenetic analysis of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] showed that this strain was a possible reassortant, as it was closely related to both porcine and porcine-like human strains, predominantly from Asia, than to typical human RVA strains. The VP6, VP7, NSP2, NSP4, and NSP5 segments of this strain showed a close similarity to porcine strains. Although the remaining gene segments (VP1, VP3, VP4, and NSP3) were closely related to human strains, all of these were porcine-like human strains [26,56,58–60,70]. With a genotype 1 (Wa-like) backbone, this finding is consistent with the hypothesis that human Wa-like strains and porcine strains have a common ancestor [5]. However, the origin of the VP2 gene of the study strain was not very definitive, as it was not only close to porcine and porcine-like human strains but also to three human strains (DC1476, DC582, and DC1127). Phylogenetically, the clusters of these three strains were shown to be distinctive from the genes of contemporary, wild-type human strains [71]. Notably, the VP7 gene of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] was located in lineage II, which comprised only porcine strains, hence implying the possibility of porcine-to-human interspecies transmission [63]. Phylogenetic analysis of porcine and human P[6] strains indicated that both porcine and human P[6] strains were present in P[6] lineages I, III, and V, hence showing that human P[6] strains might have separately emerged from at least three porcine-to-human transmissions [65]. This finding supports the Zambian G5P[6] strain, as the VP4 gene clustered and shared high nucleotide and amino acid identities with lineage V of P[6] porcine and porcine-like human strains. The NSP1 gene was most similar to porcine-like human strains. However, it was revealed to have the porcine genotype A8. Taking this together, it is likely that RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] originated by zoonotic transmission, coupled with reassortment events.

Several amino acid changes were identified in the nine variable regions when the VP7 gene of the study strain was compared to other G5 strains within each of the three lineages [64]. Additionally, the previously described conserved N-glycosylation site at residues 69–71 within the variable region 4 (VR-4) was found to be conserved in all the G5 strains used in this analysis [64,72]. Four major

antigenic regions have been described for the VP7 protein in rotaviruses (A, B, C and F) [73,74]. Marked differences in the antigenic regions of RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] were seen when it was compared to other globally circulating G5 strains. Usually, antigenic regions A and C are said to be conserved within serotypes [75]. However, multiple substitutions were observed in these regions when comparing the Zambian G5 strain to other G5 strains globally.

The amino acid sequence for the VP4 gene was 775 amino acids long and displayed amino acid identity values ranging from 91.0% to 98.3% with the reference P[6] strains. Considering it has been established that strains with amino acid identities greater than 89% belong to the same P genotype [76], our findings show that RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] belongs to the genotype P[6]. The analysis of the amino acid sequences showed that the hypervariable region (amino acid 71-208) which houses the variable region 3 (VR-3) contained most of the substitutions. Furthermore, the potential trypsin cleavage sites [68] were conserved in all the P[6] strains. Several amino acid substitutions were observed among the lineage I P[6] strains. The presence of several amino acid changes in the VP4 gene of this strain compared to other circulating P[6] strains globally is in agreement with the hypothesis that the P[6] gene has been introduced to humans via independent reassortment events [40,65,77].

Rotaviruses are genetically diverse in nature and are host-species specific, suggesting that host species barriers and restrictions exist. However, rotaviruses of animal origin may cross the host species barrier and may acquire human rotavirus gene segments, which enables the viruses to efficiently spread across human populations [4]. In this regard, G5 rotavirus strains have sporadically been documented in Latin America, Asia, Europe, and Africa [33–37,41,45,46]. Porcine P[6] strains seem to pose a lesser species barrier to humans [20]. Even though the relationship between porcine and human rotaviruses has already been established [5], whole genome analysis in this study presented the possible occurrence of interspecies transmission and reassortment between human and porcine rotaviruses.

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

#### *4.1. Ethics Statement*

This is a subset of a major project which involved the whole genome characterisation of 133 specimens collected in Zambia from 2013 through 2016 as part of the surveillance supported by the WHO/AFRO (reference 2017/757922-0) in collaboration with the University of the Free State (UFS-NGS). Ethical clearance for the main project was obtained under ethics number HSREC130/2016(UFS-HSD2016/1082) from the Health Science Research Ethics Committee (HSREC), University of the Free State, Bloemfontein, South Africa. Furthermore, this specific study was approved by the HSREC under ethics number UFS-HSD2020/0277/2104.

#### *4.2. Sample Collection*

The sample was collected in 2014 from an unvaccinated 12 month old male at Arthur Davidson Children's Hospital (ADCH) in Ndola, a rotavirus surveillance sentinel site. The child had travelled with parents from Kasama, a town in the Northern Province of Zambia which is approximately 760 km away from Ndola, Zambia. This child was admitted to a paediatric ward at ADCH, with gastroenteritis of four days duration and a history of fever. Frequency of vomiting and diarrhoea was three episodes and two episodes, respectively, in the previous 24 h. The level of dehydration was assessed as mild and the child received an oral rehydration solution and was discharged after a few days. The stool sample was screened using the enzyme immunoassay (EIA) technique for the presence of RVA antigen in the Virology laboratory in Lusaka. It was randomly picked and sent to the Diarrhoeal Pathogens Research Unit (DPRU), a World Health Organization Rotavirus Regional Reference Laboratory (WHO-RRL) in Pretoria, South Africa, as part of the WHO/AFRO annual rotavirus surveillance. Conventional genotyping was carried out at DPRU. Thereafter, the sample was shipped to the UFS-NGS unit for sequencing and whole-genome analysis.

#### *4.3. Viral dsRNA Extraction*

The viral double-stranded RNA (dsRNA) was extracted from human stool suspensions using a previously described method with modifications [78]. Approximately 100 mg stool was suspended in 200 µL phosphate-buffered saline (PBS) solution (Sigma-Aldrich®, St Louis, MO, United States). The faecal suspension was mixed with 900 µL TRI Reagent® LS (Molecular Research Centre, Cincinnati, OH, United States) and homogenized for five minutes. A 300 µL volume of chloroform (Sigma-Aldrich®, St Louis, MO, United States) was used to achieve phase separation, which was followed by centrifugation (Eppendorf microcentrifuge 5427 R, Germany) at 17,319× *g* for 20 min at 4 ◦C. The supernatant was precipitated using 700 µL ice-cold isopropanol (Sigma-Aldrich®, United States) and centrifuged (Eppendorf microcentrifuge 5427 R, Germany) at 17,319× *g* for 30 min at 4 ◦C. The supernatant was discarded, and the tubes were air-dried for 5 min, followed by the precipitation of single-stranded RNA (ssRNA) using 30 µL 8 M lithium chloride (Sigma, St Louis, MO, United States) at 4 ◦C for 16 h. The dsRNA was purified using the MinElute gel extraction kit (Qiagen, Hilden, Germany). RNA integrity was determined by electrophoresis on 1% TBE agarose gel stained with ethidium bromide (Sigma-Aldrich®, St Louis, MO, United States), which was visualised on a G: Box UV transilluminator (Syngene, Cambridge, United Kingdom).

#### *4.4. cDNA Synthesis and Purification*

cDNA synthesis was carried out using the Maxima H Minus Double-stranded cDNA kit (Thermo Fisher Scientific, Waltham, MA, United States) according to the manufacturer's instructions with minor modifications captured at the UFS-NGS SOP, whereby the dsRNA was denatured at 95 ◦C for 5 min. First strand synthesis was carried out for two hours at 50 ◦C. Random hexamer primer was employed for cDNA synthesis. The cDNA was purified using the MSB® Spin PCRapace purification kit (Stratec, Invitek Molecular, Berlin, Germany).

#### *4.5. DNA Library Preparation and Illumina*® *MiSeq Sequencing*

DNA libraries for Illumina® sequencing were prepared using the Nextera® XT DNA library preparation kit (Illumina, San Diego, CA, United States) according to the manufacturer's instructions. Briefly, DNA was tagmented at 55 ◦C for five minutes followed by ligation to Illumina® sequencing index 1 and index 2 adapters by PCR amplification. Size selection and clean-up of the DNA libraries was performed using Agencourt AMPure XP beads (Beckman Coulter, South Kraemer Boulevard Brea, CA, United States). The quantity of DNA was determined on the Qubit 2.0 fluorimeter (Invitrogen, Carlsbad, CA, United States), and a quality check of the libraries was performed on a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States). After this, sequencing was performed on an Illumina® MiSeq sequencer (Illumina, San Diego, CA, United States) using a MiSeq reagent kit v3 for 600 cycles (2 × 300 bp paired reads) with a 10% PhiX DNA control spike-in.

#### *4.6. Genome Assembly*

The raw reads obtained in FASTQ format were assembled using Geneious Prime® 2019.2.1 (https://www.geneious.com/; [79]). Briefly, the paired-end reads were merged into single reads and trimmed to remove low quality and short reads. The reads were mapped to reference sequences obtained from GenBank. Consensus sequences covering the complete open reading frame (ORF) were submitted to the National Centre for Biotechnology Information (NCBI) GenBank and assigned accession numbers MT271025–MT271035. The ORF lengths were 3267 (VP1), 2673 (VP2), 2508 (VP3), 2328 (VP4), 1194 (VP6), 981 (VP7), 1482 (NSP1), 954 (NSP2), 942 (NSP3), 528 (NSP4), and 594 (NSP5).

#### *4.7. Assignment of Genotypes*

The genotypes of each of the 11 rotavirus genome segments were determined using the online Virus Pathogen Resource (ViPR).

#### *4.8. Phylogenetic Analysis*

Gene-specific multiple sequence alignments were made using the MAFFT plugin implemented in Geneious Prime® 2019.2.1 and the MUSCLE algorithm embedded in MEGA 6.06 (for the VP2 and NSP1 segments) [80,81]. Once aligned, the DNA Model Test program in MEGA 6.06 was used to identify the optimal evolutionary model for each genome segment [82]. Using an Akaike information criterion (corrected) (AICc), the following models were found to best fit the data: HKY+G+I (VP1), GTR+G+I (VP2, VP3, and VP4), T92+G (VP6, NSP1, NSP2, NSP3, NSP4, and NSP5), and T92+G+I (VP7). Maximum likelihood trees were constructed using the optimal models in MEGA version 6.06 [82,83] with 1000 bootstrap replicates to estimate branch support [84]. The shared nucleotide and amino acid sequence identities among strains were calculated for each gene using the *p*-distance algorithm in MEGA 6.06. Analysis and visualization of the aligned concatenated whole genomes was performed on the mVISTA online platform [85].

#### **5. Conclusions**

In summary, RVA/Human-wt/ZMB/UFS-NGS-MRC-DPRU4723/2014/G5P[6] was a reassortant possessing gene segment of porcine and porcine-like human origin, and was closest to Asian strains. It is presumed that pigs play a crucial part as a source for new or newly-evolved emerging human rotaviruses. This highlights the need for continuous large-scale surveillance and whole genome analysis of circulating porcine and human rotaviruses. Furthermore, it was imperative to examine the prevalence of G5P[6] strains in Zambia. Eventually, this should result in a greater understanding of the genes that determine the transmission between hosts successfully as well as to gain insights on complex reassortment patterns between porcine and human rotaviruses.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/8/663/s1. Supplementary data 1 (S1): Identity matrices for the VP1, VP2, VP3, VP4, VP6, VP7, NSP1, NSP2, NSP3, NSP4, and NSP5 nucleotide and deduced amino acid identities among strains calculated using the p-distance algorithm in MEGA 6.06. Supplementary data 2: Comparison of amino acid sequences. Supplementary data 3: Additional phylograms of the VP1, VP2, VP3, NSP1, NSP2, NSP3, NSP4 and NSP5.

**Author Contributions:** M.M.N., M.J.M., M.L.S., and J.M.M. conceptualised the main project. W.M.M., P.N.M., J.S., E.M.M., and M.M.N. performed the laboratory experiments. M.M.N., J.S., E.M.M., M.J.M., I.P., M.L.S., and J.M.M. facilitated the obtaining of the sample. Formal analysis was performed by W.M.M., M.D.E., and M.M.N. Data curation was performed by W.M.M., P.N.M., and M.M.N. Writing (original draft preparation) was performed by W.M.M. Review of the drafts was performed by all co-authors. Supervision, project administration, and funding acquisition was conducted by M.M.N. and J.M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was principally funded by a grant awarded to M.M.N. by the World Health Organization (agreement number: UFS-AGR17-000378). Other sources of funding that supported the study were grants awarded to M.M.N.: the Bill and Melinda Gates Foundation (BMGF-OPP1180423\_2017), South African Medical Research Foundation through the Self-Initiated Research grant (SAMRC-SIR), National Research Foundation (NRF-120814), and a grant from the Poliomyelitis Research Foundation (PRF-19/95).

**Acknowledgments:** We greatly thank the Zambia team who assisted with sample collection and ELISA testing. We would like to thank Khutso Mothapo, Kebareng Rakau, and Nonkululeko Magagula at the WHO-RRL (Pretoria, South Africa) who assisted with the retrieval of samples. We would also like to thank Sebotsana Rasebotsa, Milton Mogotsi, Emmanuel Ayodeji, Teboho Mooko, and Gilmore Pambuka for their assistance with laboratory work. We are grateful to Stephanus Riekert for technical ICT support.

**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 World Health Organization.

#### **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*

### **Phylogenetic Analyses of Rotavirus A from Cattle in Uruguay Reveal the Circulation of Common and Uncommon Genotypes and Suggest Interspecies Transmission**

**Matías Castells 1,2,\* , Rubén Darío Ca**ff**arena 2,3 , María Laura Casaux <sup>2</sup> , Carlos Schild <sup>2</sup> , Samuel Miño <sup>4</sup> , Felipe Castells <sup>5</sup> , Daniel Castells <sup>6</sup> , Matías Victoria <sup>1</sup> , Franklin Riet-Correa <sup>2</sup> , Federico Giannitti <sup>2</sup> , Viviana Parreño <sup>4</sup> and Rodney Colina 1,\***


Received: 7 April 2020; Accepted: 30 June 2020; Published: 14 July 2020

**Abstract:** Uruguay is one of the main exporters of beef and dairy products, and cattle production is one of the main economic sectors in this country. Rotavirus A (RVA) is the main pathogen associated with neonatal calf diarrhea (NCD), a syndrome that leads to significant economic losses to the livestock industry. The aims of this study are to determine the frequency of RVA infections, and to analyze the genetic diversity of RVA strains in calves in Uruguay. A total of 833 samples from dairy and beef calves were analyzed through RT-qPCR and sequencing. RVA was detected in 57.0% of the samples. The frequency of detection was significantly higher in dairy (59.5%) than beef (28.4%) calves (*p* < 0.001), while it did not differ significantly among calves born in herds that were vaccinated (64.0%) or not vaccinated (66.7%) against NCD. The frequency of RVA detection and the viral load were significantly higher in samples from diarrheic (72.1%, 7.99 log<sup>10</sup> genome copies/mL of feces) than non-diarrheic (59.9%, 7.35 log<sup>10</sup> genome copies/mL of feces) calves (*p* < 0.005 and *p* = 0.007, respectively). The observed G-types (VP7) were G6 (77.6%), G10 (20.7%), and G24 (1.7%), while the P-types were P[5] (28.4%), P[11] (70.7%), and P[33] (0.9%). The G-type and P-type combinations were G6P[11] (40.4%), G6P[5] (38.6%), G10P[11] (19.3%), and the uncommon genotype G24P[33] (1.8%). VP6 and NSP1-5 genotyping were performed to better characterize some strains. The phylogenetic analyses suggested interspecies transmission, including transmission between animals and humans.

**Keywords:** rotavirus; bovine; genotypes; interspecies transmission; diarrhea

#### **1. Introduction**

Neonatal calf diarrhea (NCD) is a syndrome of worldwide distribution and the major cause of mortality of dairy calves before weaning [1]. NCD has a negative impact on animal welfare and leads to significant economic losses to the livestock industry [2–5].

Rotavirus A (RVA) is the main pathogen associated with NCD [6,7]. RVA (species *Rotavirus A*; genus *Rotavirus*; subfamily *Sedoreovirinae*; family *Reoviridae*) is a nonenveloped virus with a triple-layered capsid and a genome composed of 11 segments of double-stranded RNA [8]. RVA is widespread in dairy farms in Uruguay, and viable viral particles have been detected in sources of drinking water used for calves [9], suggesting water contamination and waterborne transmission.

Rotaviruses are classified by a binary system of G and P types for VP7 and VP4, respectively, determined by sequence analyses. In 2008, a complete genome classification system, named genotype constellation, assigning a specific genotype to each of the 11 genome segments was developed [10]. The VP7-VP4-VP6-VP1-VP2-VP3-NSP1-NSP2-NSP3-NSP4-NSP5/6 genes of rotavirus strains are classified using the abbreviations Gx-P[x]-Ix-Rx-Cx-Mx-Ax-Nx-Tx-Ex-Hx (where x is the genotype number), respectively.

Recently, since the inclusion of gene segments other than VP7 and VP4 in molecular analyses, gene reassortment has been described as a common event in RVA, sometimes between virus strains originated from different hosts, suggesting interspecies transmission [10–13].

Surveys describing the epidemiology of RVA in cattle in South America are mainly restricted to Brazil and Argentina; no published data about RVA epidemiology in Uruguayan calves are available. However, other viruses such as bovine coronavirus and bovine astrovirus have been detected in Uruguay [14,15].

Uruguay is one of the main exporters of beef [16] and dairy products [17]. Furthermore, cattle production is one of the main economic sectors in this country, with almost 12 million head of cattle accounting for 33% of the total exports [18]. The aims of this study are to determine the frequency of RVA infections and to analyze the genetic diversity of the RVA strains detected in Uruguayan calves.

#### **2. Results**

#### *2.1. Detection Frequency of RVA in Uruguayan Calves*

Rotavirus A was detected in 57.0% (475/833) of the analyzed samples. The frequency of detection was significantly higher in dairy (59.5%, 456/766) than beef (28.4%, 19/67) calves (OR: 3.72, 95% CI: 2.14–6.44; *p* < 0.000001; Figure 1a). The frequency of RVA detection in live calves was higher (58.0%, 444/766) than in deceased calves (46.3%, 31/67), although this difference was not statistically significant (*p* = 0.06; Figure 1b). The frequency of detection in dairy calves born in herds that vaccinated (64.0%, 144/225) or did not vaccinate dams (66.7%, 164/246) against NCD did not differ significantly (*p* = 0.5; Figure 1c). The frequency of RVA detection was significantly higher in samples from diarrheic (72.1%, 173/240) than non-diarrheic (59.9%, 163/272) dairy calves (OR: 1.73, 95% CI: 1.19–2.50; *p* < 0.005; Figure 1d). No seasonal distribution was observed in RVA detection (data not shown).

Rotavirus A was detected in 58.8% (87/148), 70.6% (142/201), 68.2% (75/110), and 52.9% (18/34) of dairy calves in the first, second, third, and fourth weeks of life, respectively (Table 1). Statistically significant differences were observed between the second and the first weeks of age (OR: 1.69, 95% CI: 1.08–2.64; *p* = 0.02), and between the second and the fourth weeks of age (OR: 2.14, 95% CI: 1.02–4.48; *p* = 0.04). The mean age in days of RVA-positive dairy calves was significantly lower in diarrheic than nondiarrheic calves (*p* = 0.02; Table 1).

The RVA viral load was significantly higher in diarrheic than nondiarrheic dairy calves (*p* = 0.007; Table 1), ranging between 1.14 × 10<sup>4</sup> and 7.36 × 10<sup>12</sup> genome copies/milliliter (gc/mL) of feces. In all four age groups, the frequency of RVA detection was higher in diarrheic than nondiarrheic dairy calves: 69.0% (40/58) vs. 52.2% (47/90) in the first week, 72.1% (98/136) vs. 67.7% (44/65) in the second week, 68.8% (22/32) vs. 67.9% (53/78) in the third week, and 85.7% (6/7) vs. 44.4% (12/27) in the fourth week

of age. A statistically significant difference was observed only within the first week (OR: 2.03, 95% CI: 1.01–4.07; *p* = 0.04).

**Figure 1.** Frequency of Rotavirus A (RVA) detection in calves. (**a**) Frequency of RVA detection in dairy vs. beef calves; (**b**) frequency of RVA detection in live vs. deceased calves; (**c**) frequency of RVA detection in calves from vaccinated <sup>a</sup> vs. unvaccinated dairy herds; (**d**) frequency of RVA detection in diarrheic vs. non diarrheic dairy calves. Comparisons with statistically significant differences are indicated. <sup>a</sup> Most of the vaccines against neonatal calf diarrhea available in Uruguay include two RVA strains.



<sup>a</sup> Mean age in days of RVA-positive calves. <sup>b</sup> Mean RVA viral load expressed as log10 of RVA genome copies per milliliter of feces. Equal numbers in superscript refer to values with statistically significant differences (*p* < 0.05).

#### *2.2. VP7 and VP4 Genotyping*

We obtained 58 and 116 sequences for VP7 and VP4, respectively. The detected G-types (VP7) were G6 (77.6%, 45/58), G10 (20.7%, 12/58), and G24 (1.7%, 1/58), while the P-types (VP4) were P[5] (28.4%, 33/116), P[11] (70.7%, 82/116), and P[33] (0.9%, 1/116). The following G- and P-type combinations were obtained for 57 strains: G6P[11] (40.4%, 23/57), G6P[5] (38.6%, 22/57), G10P[11] (19.3%, 11/57), and G24P[33] (1.8%, 1/57). Furthermore, 60 strains had undetermined G- or P-type: GXP[11] (80.0%, 48/60), GXP[5] (18.3%, 11/60), and G10P[X] (1.7%, 1/60).

#### *2.3. VP6 and NSP1-5 Genotyping*

Ten samples, including representative VP7 and VP4 genotype combinations observed, were selected for VP6 and NSP1-5 gene characterization: 2 G6P[5], 2 G6P[11], 2 G10P[11], 2 GXP[11], 1 G10P[X], and 1 G24P[33] (Table 2). All the strains were I2 (VP6), N2 (NSP2), and E12 (NSP4). Nine were H3 and one could not be determined HX (NSP5). Five strains were A3, four were A13, and one could not be determined AX (NSP1). Eight strains were T6, one was T9, and one could not be determined TX (NSP3).


**Table 2.** Genotype constellation of 10 RVA strains from Uruguayan calves.

Uncommon genotypes are shadowed in grey.

#### *2.4. Phylogenetic Analyses*

The phylogenetic analyses showed an intricate genetic scenario. The analyses of the VP7 gene showed that G6 and G10 Uruguayan strains clustered in two and one different lineages, respectively, with sequences obtained from cattle. Specifically, the G6P[5] Uruguayan strains clustered in one lineage (split into two sublineages) with Argentinian strains, and the G6P[11] Uruguayan strains clustered separately in a lineage with Slovenian strains (Figure 2). The G10 Uruguayan strains clustered in a lineage (split into two sublineages) with Argentinian strains (Figure 3). Brazilian G6 and G10 strains clustered separately with Uruguayan and Argentinian G6 and G10 strains.

The phylogenetic analyses of the VP4 gene showed that P[5] Uruguayan strains clustered in a lineage with Argentinian G6P[5] strains obtained from cattle, and Brazilian P[5] strains clustered separate (Figure 4). The P[11] Uruguayan strains clustered in three lineages with sequences obtained from cattle, two of the lineages were comprised of G6 and G10 Argentinian strains (and one of these lineages is split into two sublineages), and the other lineage comprised of G6P[11] Brazilian strains, although P[11] Uruguayan strains were distinct to the majority of the Brazilian P[11] strains (Figure 5).

In the phylogenetic tree of the NSP1 gene, we observed that Uruguayan strains clustered in three different genetic lineages of the genotype A3: one jointly with human strains from Paraguay and Brazil, another with Italian and Belgian human strains, and another with a goat strain from Argentina and, in one genetic lineage of the genotype A13, with an Argentinian strain from a cow (Figure S1).

The phylogenetic analysis of the NSP2 gene showed that the Uruguayan strains were clustered in two separate lineages: one with Argentinian strains from cow and goat, and the other with strains from guanaco and vicuña from Argentina and strains from humans from Australia (Figure S2).

On the other hand, the phylogenetic analysis of the NSP3 gene showed that the T6 Uruguayan strains were clustered in three sublineages within one lineage: one together with strains distributed worldwide (including vaccine strains), one with Argentinian (vicuña and guanaco), Japanese (cow), Slovenian (human), and Paraguayan (human) strains, and the third with a goat strain from Argentina and a human strain from Belgium. The T9 strain clustered with the other four T9 strains detected so far (from Japan and the USA; Figure S3).

For the NSP4 gene, we observed that besides the Uruguayan strains obtained in our study, only sequences from South America were available. The phylogenetic analysis showed that Uruguayan strains clustered in four different lineages together with strains from several host species (cows, guanacos, horses, goats, and humans), all from this subcontinent (Figure S4).

The phylogenetic analysis of the NSP5 gene showed that the Uruguayan strains were clustered in three sublineages within one lineage: one together with strains distributed worldwide in several host species), other with an Argentinian strain from a cow and a Paraguayan strain obtained from a human, and another with a strain from a guanaco from Argentina, a strain from a yak from China, and a strain from a human from Hungary (Figure S5).

Lastly, the phylogenetic analysis of the VP6 gene showed that the Uruguayan strains were clustered in three lineages: one conformed only with Uruguayan strains, another lineage with an Argentinian strain from a cow, and another lineage with South American strains from various hosts (human, llama, sheep, and goat), Japanese strains from human and cow, and a roe deer Slovenian strain (Figure S6).

**Figure 2.** Maximum likelihood tree of the G6 genotype of the VP7 gene. The best nucleotide substitution model (TIM2 + I + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in red. Shimodaira–Hasegawa-approximate likelihood-ratio test (SH-aLRT) values ≥ 80 are shown.

**Figure 3.** Maximum likelihood tree of the G10 genotype of the VP7 gene. The best nucleotide substitution model (TPM3 + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in red. SH-aLRT values ≥ 80 are shown.

**Figure 4.** Maximum likelihood tree of the P[5] genotype of the VP4 gene. The best nucleotide substitution model (TIM + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in red. SH-aLRT values ≥ 80 are shown.

**Figure 5.** Maximum likelihood tree of the P[11] genotype of the VP4 gene. The best nucleotide substitution model (TPM3u + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in red. SH-aLRT values ≥ 80 are shown.

#### **3. Discussion**

Rotavirus A was detected in feces and intestinal contents collected from dairy and beef calves with a frequency of 57%, which was higher than reports from Argentina and Brazil (17–42%) [19–22], and other geographic regions (20–49%) [7,23–25]. On the other hand, in Australia, the frequency of RVA detection was 80%, which is higher than the detected in our study [6]. Interestingly, most of the mentioned studies were conducted by assays different than RT-qPCR, except the one conducted in Australia. It is well documented that the RT-qPCR for RVA detection has a higher sensitivity than other assays, reducing the risk of false-negatives (i.e., ELISA, electron microscopy, PAGE, immunochromatography, and conventional PCR) [6,26–28], which could explain the higher frequency observed in Uruguay when compared with neighboring countries while reducing the risk of false-positive results, also given its higher specificity. Furthermore, the use of RT-qPCR, which is known to detect very few genomic copies, allows pathogen detection in clinical and subclinical calves. In addition, in many field situations, the time of onset of diarrhea is not known, so the peak of pathogen shedding may have already passed, or the infection could be just settling down by the time of sampling [29]. The limit of detection in our study (10<sup>4</sup> gc/mL of feces) and the higher RVA viral load in diarrheic than nondiarrheic calves are in agreement with the stated by Torres-Medina et al. [29]. On the other hand, we also observed high viral loads in some nondiarrheic calves.

Infection with RVA has long been associated with diarrhea [29–31], as observed in our study, where RVA detection was more frequent in diarrheic than in nondiarrheic calves, independently of their age (up to 4 weeks). Concerning the calves' age, we observed that the proportion of calves shedding RVA was higher in the second and third weeks of age, as observed in Brazil [19,32] and elsewhere [33]. In addition, the mean age of RVA-positive calves in our study is similar to the age reported previously [31], and we observed that diarrheic calves positive for RVA were younger than nondiarrheic calves, indicating that calves are exposed to this pathogen early after birth.

Although the sampling between beef and dairy farms was unequal, our results indicate that the circulation of RVA was higher in dairy than beef calves. This contrasts with the reported results in neighboring countries, where RVA was more frequently detected in beef than dairy calves [19,20] or in a similar frequency [21]. Our results also contrast with those observed in a study conducted in Australia [6].

A common practice used to prevent NCD is the vaccination of pregnant cows/heifers during the last stage of pregnancy to protect the calves by the transference of passive maternal antibodies through colostrum intake. Most available vaccines in the Uruguayan market include bovine rotavirus A strains (most of them include two strains, G6 and G10, as detailed by the manufacturers). In this study, we observed a similar frequency of RVA detection in calves from vaccinated and unvaccinated herds. Failure in the protection against RVA infection by the vaccine was reported in studies conducted in Argentina and Brazil [34–37]; although vaccines are not effective in preventing RVA infection, they significantly reduce morbidity, the severity of diarrhea, and mortality related to RVA [38].

In this study, we determined the RVA genotypes circulating in calves in Uruguay. Overall, the VP7 and VP4 genotypes observed in this country are the most prevalent in cattle worldwide [39], although, unexpectedly, we detected a G24P[33] strain, which thus far had only been reported from an asymptomatic cow and her calf in Japan [11]. The G24P[33] strain detected in Uruguay was obtained from a 10-day-old asymptomatic dairy calf sampled in August 2016.

Regarding the VP6 and NSP1-5 genotyping, the Uruguayan strains, including the G24P[33], showed a relatively conserved genotype constellation I2-A3/A13-N2-T6/T9-E12-H3, corresponding to VP6 and NSP1-5 genotypes, respectively. These genotypes are commonly found in cattle, with the exception of T9 [40]. The T9 genotype has been sporadically detected in two cows from Japan [11], in a child from Japan [41], and in a child from the USA [42]. This genotype has been associated with atypical VP7 and VP4 genotypes (G21P[29], G24P[33], G8P[14], and G24P[14]). In this study, we observed the T9 genotype associated with G24P[33]. Indepth analysis of the RVA/Cow-wt/URY/LVMS3024/2016/G24P[33] strain revealed almost the same genotype constellation as the RVA/Cow-wt/JPN/Dai-10/2007/G24P[33] strain from Japan, with the unusual G24, P[33], and T9 genotypes. The only difference was observed in the NSP4 gene that was E12 in the Uruguayan strain and E2 in the Japanese. It is interesting to note that all the Uruguayan strains were E12, a genotype widely detected in cattle [12], guanacos [12], horses [43,44], goats [45], and children [46,47] in South America. This reinforces the notion that the E12 genotype may be restricted to South America, as previously postulated [44].

The rare G24P[33] strain detected in our study represented a challenge. The G24, P[33], and T9 genotypes observed in this strain provides information for a possible introduction of the virus from Japan to Uruguay, or vice versa. The expansion of the Wagyu beef industry beyond Japan [48] could have influenced the dispersion of some RVA strains through live cattle exports. On the other hand, the E12 genotype in the Uruguayan G24P[33] strain and E2 genotype in the Japanese G24P[33] strain represented a probable gene reassortment, which is a more plausible scenario than the emergence of two independent strains with the same rare genotype constellation except for NSP4. Further studies should be conducted to determine the evolution and possible emergence of these rare genotypes.

In the phylogenetic analyses of all the genes, it can be observed that Uruguayan strains clustered mainly with South American strains. The only gene that did not show any South American-specific lineage was NSP3, in which the Uruguayan strains clustered mainly with Argentinian strains, but also with strains from other continents. These data, together with the identification of the E12 genotype in all the Uruguayan sequences, suggest a South American origin of RVA lineages [44]. Furthermore, the phylogenetic analyses showed an intricate pattern of diversity, with evidence of gene reassortments, interspecies transmission, local dispersion of some strains, and circulation of strains that are most prevalent in cattle worldwide.

The analyses of VP7 and VP4 showed a conserved pattern with all the Uruguayan strains clustering, with strains detected only in cattle and mainly from Argentina, indicating a probable host species and geographic linkage. Due to the shortage of G24 and P[33] sequences in the database (2 and 1, respectively), no phylogenetic analyses were performed for these genotypes. In the VP7 and VP4 phylogenetic analysis, the majority of strains characterized in this study clustered closely with strains detected in Argentinian cattle. The exceptions were one G6 lineage that clustered with European strains isolated from cattle, and one in P[11] sublineage that clustered with Brazilian strains isolated from cattle. There is a clear phylogenetic relationship between the strains detected in the cattle in Uruguay and Argentina, whereas Brazilian strains were, in general, phylogenetically distant from the Uruguayan strains. In addition, Uruguayan strains clustered together among themselves, suggesting that limited introductions of RVA into the country have occurred, but the strains were widely dispersed in the cattle. A possible explanation for the genetic similarity between the Uruguayan and Argentinian strains and their divergence to the Brazilian strains could be explained, in part, by the breed of cattle. In Uruguay and Argentina, most of the cattle breeds are *Bos taurus*, while in Brazil, there are mostly *Bos indicus* or *Bos indicus* x *Bos taurus* crosses. Although it has not been studied in cattle, different human subpopulations appeared to have different susceptibility infection and clinical disease, and this susceptibility is dependent on the rotavirus genotype, and in some cases, it also depends on different rotavirus strains of the same genotype [49].

Based on the phylogenetic analyses, we observed evidence of gene reassortment and interspecies transmission events. Regarding the former event, in addition to the previously mentioned gene reassortment of the G24P[33] strain, strong evidence was observed in the strains RVA/Cow-wt/URY/ LVMS1812/2016/G6P[5] and RVA/Cow-wt/URY/LVMS3206/2016/GxP[11] because both strains clustered together in all the genes, except in VP4 (which showed different genotypes, (P[5] and P[11], respectively), indicating that a possible gene reassortment event may have occurred. Another piece of evidence was observed in the RVA/Cow-wt/URY/LVMS1788/2016/GxP[11] strain because it clustered together with other Uruguayan strains in most of the genes, except in NSP1 and NSP3 genes, which clustered alone in different genetic lineages, also suggesting a gene reassortment event. Furthermore, the strain RVA/Cow-wt/URY/LVMS1837/2016/G10P[11] clustered together with RVA/Cow-wt/URY/LVMS2625/ 2016/G10P[11] and RVA/Cow-wt/URY/LVMS3053/2016/G10P[x] in most of the genes, but clustered

separately in distant genetic lineages in NSP2 and NSP5; this was probably due to gene reassortment. On the other hand, an interesting observation was that, in general, G6 strains tended to cluster together in most of the genes, and the same was observed for the G10 strains, with the exceptions aforementioned.

Regarding interspecies transmission, we observed that in the analyses of VP7 and VP4, all the Uruguayan strains clustered with other bovine strains, so these gene segments seem to be more host-specific than the other genes. On the other hand, and based on the phylogenetic analyses, we observed evidence suggesting interspecies transmission because the bovine strains detected in Uruguay closely clustered with strains detected in other host species. We observed that bovine Uruguayan strains A13 (NSP1 gene) clustered together with strains isolated from humans and a goat, possibly indicating events of interspecies transmission. Two lineages showed a close relationship between Uruguayan bovine strains and human strains (from South America and Europe); these human strains were reported to be Artiodactyl-like and a product of interspecies transmission [10,47,50], as well as the goat strain of a third lineage [45], which is in accordance with our results. In the NSP2-5 and VP6 genes, we observed that the Uruguayan bovine strains clustered in some lineages with strains isolated from other host species (human, goat, guanaco, vicuna, roe deer, llama, and sheep), mainly from South America, that were proposed to be originated by interspecies transmission [12,45,47,51], again in accordance with our results. Another piece of evidence supporting this event was observed in the NSP4; all the RVA strains detected in South America were E12, independent of the host species where they were isolated (horse, cow, guanaco, human, goat), suggesting interspecies transmission and fixation of this genotype in South America [44]. The interspecies transmission of RVA is widely documented [10–13], and our results support this event. In South America, it is common to raise different livestock species on the same farm in close contact with humans [45], which increases the possibility of interspecies transmission. Our results support that interspecies transmission is a common event in South America, including the possibility of zoonotic transmission [45,51,52].

Lastly, our study had some limitations. In Uruguay, dairy farming is concentrated in the southwest region and calves are raised under intensive production systems that facilitated the collection of the samples, while beef calves are mostly bred in extensive production systems and dispersed throughout the country, which hindered the access to samples. This resulted in an overrepresentation of dairy (92%) versus beef (8%) samples in our study. Another limitation was that we had no spiked control to determine if there was inhibition of the qPCR, which may lead to false-negatives. Regarding coinfections, the methodology used has the limitation that sequences obtained from a single animal would have only represented the predominant strain and/or sequences with multiple traces that were not included in the study. It is important to mention that, from our analyses, we could not determine the route nor the time in which the gene reassortment and the interspecies transmission events took place.

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

#### *4.1. Samples*

Fecal samples of 766 live calves and intestinal contents from 67 naturally-deceased calves were collected from 833 different calves from dairy and beef herds in Uruguay between 2015 and 2018. Sampled herds were distributed in 10 of the 19 regions of the country (Figure 6), and throughout the year, including samples collected in the four climate seasons. In addition, 766 samples were from dairy calves, and 67 from beef calves. We compared the frequency of the RVA infection between groups only for dairy calves. A total of 240 dairy calves had diarrhea at the time of sampling, while 272 were nondiarrheic dairy calves (this information was unavailable for 321 calves). The distribution by age in the first, second, third, and fourth weeks of life was 148, 201, 110, and 34 dairy calves, respectively (the age was unavailable for 340 calves). A total of 225 calves were from dairy herds vaccinated against NCD and 246 calves were from nonvaccinated dairy herds (herd vaccination history was unavailable for 362 calves).

**Figure 6.** Map of Uruguay, the regions from which samples were collected shown in grey.

#### *4.2. Sample Suspension, RNA Extraction, Reverse Transcription, Detection and Quantification of RVA*

Samples were diluted 1:10 (*v*:*v*) in phosphate-buffered saline solution, centrifuged at 3000× *g* for 20 min at 4 ◦C, and supernatants were collected and stored at −80 ◦C. Viral RNA was extracted using a QIAamp® cador® Pathogen Mini Kit (Qiagen®, Hilden, Germany), following the manufacturer's instructions. Reverse transcription (RT) was carried out with RevertAid® Reverse Transcriptase (Thermo Fisher Scientific®, Waltham, MA, USA) and random hexamers primers (Qiagen®), following the manufacturer´s instructions. All RNAs and cDNAs were stored at −80 ◦C until further viral analyses. Screening and quantification of the samples for RVA identification were carried out through a quantitative polymerase chain reaction (qPCR) targeted to the NSP3 gene, as described elsewhere [9]. Briefly, 12.5 µL of SensiFAST™ Probe No-ROX Kit (Bioline®, London, UK), 5.0 µL of nuclease-free water, 1.0 µL of 10 µM forward primer, 1.0 µL of 10 µM reverse primer, 0.5 µL of 10 µM probe, and 5 µL of cDNA were mixed in 0.2-mL PCR tubes. All samples were analyzed in duplicate. In order to validate the complete process, an RVA-positive (G6P[5] strain) and an RVA-negative fecal sample were used as positive and negative controls, respectively.

#### *4.3. Rotavirus A Genotyping*

Quantitative-PCR positive samples were subsequently subjected to amplification of VP7 and VP4 (VP8\*). Briefly, 12.5 µL of MangoMix™ (Bioline®), 5 µL of cDNA, 5.5 µL of nuclease-free water, 1 µL of dimethyl sulfoxide, 0.5 µL of 20 µM forward primer and 0.5 µL of 20 µM reverse primer were mixed in 0.2-mL PCR tubes. Forward and reverse primers for VP7 and VP4 (VP8\*) amplification are described elsewhere [20,53]. In addition, 10 samples, including representative VP7 and VP4 genotype combinations observed in this study, were selected for VP6 and NSP1-5 gene characterization. Primers and cycling conditions were used, as described elsewhere [10], and PCR reagents were used, as described above. Genotyping was performed using the web-based genotyping tool RotaC v2.0 [54].

#### *4.4. PCR Product Purification, Sequencing, and GenBank Accession Numbers*

PCR products were visualized in 1–2% agarose gels and positive samples were purified using PureLink™ Quick Gel Extraction and PCR Purification Combo Kit (Invitrogen®, Carlsbad, CA, USA), according to the manufacturer's instructions. Both cDNA strands were sequenced by Macrogen Inc. (Seoul, Korea). Sequences were deposited in GenBank with accession numbers: MN649559—MN649674 (VP4), MN649675—MN649732 (VP7), MN649733—MN649742 (VP6), MN649743—MN649751 (NSP1), MN649752—MN649761 (NSP2), MN649762—MN649770 (NSP3), MN649771—MN649780 (NSP4), and MN649781—MN649789 (NSP5).

#### *4.5. Phylogenetic Analysis*

All the available sequences corresponding to the genotypes observed in the RVA strains detected in this study, previously determined with RotaC, were downloaded from the Virus Variation Resource (http://www.ncbi.nlm.nih.gov/genome/viruses/variation/) [55]. A dataset was created for each genotype, and multiple sequence alignments were obtained using ClustalW implemented in MEGA 7 software [56]. The final alignment of each gene comprised all the worldwide sequences that covered the length of the sequences obtained in this study. The length of the sequences and the nucleotide position, involved in the phylogenetic analysis of each gene, are detailed in Table 3. The nucleotide substitution models that best fit each dataset (Table 3) and the maximum likelihood trees were obtained using W-IQ-TREE (available at http://iqtree.cibiv.univie.ac.at) [57]. The branches support was estimated with the Shimodaira–Hasegawa-approximate likelihood-ratio test (SH-aLRT) [58]. Trees were visualized in FigTree (http://tree.bio.ed.ac.uk/software/figtree/).

**Table 3.** Information about the final alignments obtained for the phylogenetic analyses.


\* Reference strain: WC3.

#### *4.6. Statistical Analyses*

Data were organized and graphics were generated using Microsoft® Office Excel. Categorical data were evaluated with RStudio v1.0.136 software through Pearson's chi-squared tests. Odds ratios (OR) and 95% confidence intervals (CI) were calculated with jamovi software (available at https://www. jamovi.org/). Viral load values (genome copies/milliliter of feces) were log10 transformed. For the viral load and mean age analyses, the Shapiro–Wilk test was performed, rejecting the normality of the data, so the Mann–Whitney U test was performed with the same software. For all tests, differences were considered statistically significant if the obtained *p*-value was < 0.05.

#### **5. Conclusions**

Rotavirus A is widespread in cattle in Uruguay and is associated with diarrhea in calves, with a peak of viral shedding at 2–3 weeks of age, and higher viral shedding in diarrheic versus non-diarrheic calves. Even though the main genotypes observed in this country are the most prevalent worldwide, a rare strain was detected with a G24-P[33]-I2-A13-N2-T9-E12-H3 genotype constellation. The E12 genotype detected in all strains, regardless of the VP7 and VP4 genotypes, appears to be a South American geographic marker. An intricate genetic scenario was evidenced, with gene reassortment and interspecies transmission events, including transmission between animals and humans.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/7/570/s1, Figure S1: Maximum likelihood tree of the NSP1 gene. The best nucleotide substitution model (TIM + I + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in different colors. SH-aLRT values ≥ 80 are shown. Figure S2: Maximum likelihood tree of the NSP2 gene. The best nucleotide substitution model (TIM + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in different colors. SH-aLRT values ≥ 80 are shown. Figure S3: Maximum likelihood tree of the NSP3 gene. The best nucleotide substitution model (TIM3 + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in different colors. SH-aLRT values ≥ 80 are shown. Figure S4: Maximum likelihood tree of the NSP4 gene. The best nucleotide substitution model (HKY + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in different colors. SH-aLRT values ≥ 80 are shown. Figure S5: Maximum likelihood tree of the NSP5 gene. The best nucleotide substitution model (TN + I + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in different colors. SH-aLRT values ≥ 80 are shown. Figure S6: Maximum likelihood tree of the VP6 gene. The best nucleotide substitution model (TIM + I + G) and the maximum likelihood tree were obtained with W-IQ-TREE. Uruguayan strains are shown in different colors. SH-aLRT values ≥ 80 are shown.

**Author Contributions:** Conceptualization, M.C. and R.C.; methodology, M.C., R.D.C., M.L.C., C.S., S.M., F.C. and D.C.; resources, M.C., F.R.-C., F.G., V.P. and R.C.; writing—original draft preparation, M.C.; writing—review and editing, M.C., R.D.C., M.L.C., S.M., F.C., D.C., M.V., F.R.-C., F.G., V.P. and R.C.; funding acquisition, M.C., F.R.-C. and R.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by "Instituto Nacional de Investigación Agropecuaria" (INIA), grant number PL\_015 N-15156 and N-23398, and by the "Universidad de la República" program "Polo de Desarrollo Universitario". The APC was funded by "Universidad de la República".

**Acknowledgments:** M.C. acknowledges support from the "Agencia Nacional de Investigación e Innovación" (ANII) through a PhD scholarship, and "Comisión Sectorial de Investigación Científica" (CSIC) and ANII for mobility fellowships.

**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.

#### **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* **Group A Rotavirus Detection and Genotype Distribution before and after Introduction of a National Immunisation Programme in Ireland: 2015–2019**

#### **Zoe Yandle \* , Suzie Coughlan, Jonathan Dean, Gráinne Tuite, Anne Conroy and Cillian F. De Gascun**

UCD National Virus Reference Laboratory, University College Dublin, Dublin 4, Ireland; suzie.coughlan@ucd.ie (S.C.); jonathan.dean@ucd.ie (J.D.); grainne.tuite@ucd.ie (G.T.); anne.conroy@ucd.ie (A.C.); cillian.degascun@ucd.ie (C.F.D.G.) **\*** Correspondence: zyandle@ucd.ie; Tel.: +353-1-716-4401

Received: 19 May 2020; Accepted: 5 June 2020; Published: 7 June 2020

**Abstract:** Immunisation against rotavirus infection was introduced into Ireland in December 2016. We report on the viruses causing gastroenteritis before (2015–2016) and after (2017–2019) implementation of the Rotarix vaccine, as well as changes in the diversity of circulating rotavirus genotypes. Samples from patients aged ≤ 5 years (n = 11,800) were received at the National Virus Reference Laboratory, Dublin, and tested by real-time RT-PCR for rotavirus, Rotarix, norovirus, sapovirus, astrovirus, and enteric adenovirus. Rotavirus genotyping was performed either by multiplex or hemi-nested RT-PCR, and a subset was characterised by sequence analysis. Rotavirus detection decreased by 91% in children aged 0–12 months between 2015/16 and 2018/19. Rotarix was detected in 10% of those eligible for the vaccine and was not found in those aged >7 months. Rotavirus typically peaks in March–May, but following vaccination, the seasonality became less defined. In 2015–16, G1P[8] was the most common genotype circulating; however, in 2019 G2P[4] was detected more often. Following the introduction of Rotarix, a reduction in numbers of rotavirus infections occurred, coinciding with an increase in genotype diversity, along with the first recorded detection of an equine-like G3 strain in Ireland.

**Keywords:** gastroenteritis; rotavirus; Rotarix; pediatric; diagnostics; molecular epidemiology; G3P[8]; equine-like

#### **1. Introduction**

Rotavirus is a leading cause of pediatric acute gastroenteritis, causing fever, vomiting, and diarrhoea. Mortality rates are highest in low income developing countries, where it causes approximately 128,000 fatal cases per year in those under five years old [1,2]. With the availability of rotavirus vaccines, the rate of global hospitalisations due to rotavirus or acute gastroenteritis, as well as deaths due to acute gastroenteritis, has decreased [3]. In Europe, pediatric rotavirus infection results in approximately 75,000−150,000 hospitalisations annually, with 2−4 times more children seeking out-patient medical care [4]. In Ireland, average crude incidence rates were 55 per 100,000 population in the 2007–2015 period [5], with hospitalisation rates of approximately 1190 per 100,000 [6], compared to the majority of EU member states, which report rates of 300–600 per 100,000 [4].

In 2009, the World Health Organisation (WHO) recommended global rotavirus vaccination [7] and, in Europe, 13 countries include it in their universal immunisation programmes, with a further five offering the vaccine for certain risk groups, specific regions, or requiring partial payment [8]. Two licensed live-attenuated vaccines are available in Europe; the pentavalent bovine-human reassortment rotavirus vaccine, RotaTeq (Merck & Co., West Point, PA, USA), and the human monovalent vaccine, Rotarix (GlaxoSmithKline, Rixensart, Belgium). In December 2016, Rotarix was introduced into the Irish national immunisation programme, with the vaccine administered in two doses at 2 and 4 months of age. Most recent figures (Q3, 2019) show the national uptake of the vaccine is 89% [9]. Rotavirus is a notifiable disease in Ireland, and laboratory confirmed cases are reported to the Health Protection Surveillance Centre. Effectiveness of both RotaTeq and Rotarix has been well documented, with the UK, Germany, and Belgium reporting an approximate 85% reduction in the presentation of severe rotavirus disease following vaccination [10–12].

Rotaviruses are double stranded RNA viruses containing 11 genome segments. There are 10 groups, A–J, defined by the middle VP6 capsid antigen, [13] two of which (I and J) were recently discovered in dogs and bats, respectively [14,15]. However, in humans, the majority of infections are caused by Group A rotavirus. Classification is a binary system depending on the expression of two outer proteins; the G and P-type, encoded by VP7 and VP4, respectively. Full genome analysis (where the VP7-VP4-VP6-VP1-VP2-VP3-NSP1-NSP2-NSP3-NSP4-NSP5/6 genes of rotavirus (RV) strains are described using the abbreviations Gx-P[x]-Ix-Rx-Cx-Mx-Ax-Nx-Tx-Ex-Hx), is required to monitor the evolution of the virus and detect reassortment [16].

Despite the theoretical possibility for numerous rotavirus G/P constellations, six account for 80–90% of circulating genotypes, namely G1P[8], G2P[4], G3P[8], G4P[8], G9P[8] and G12P[8]. Distribution of these commonly detected genotypes can vary by year, country, and age [17]. Despite the natural fluctuation of genotype diversity, increasing data suggest that the changes may be due to the impact of strain-specific vaccines [18]. Both in Belgium and the UK, before immunisation, G1P[8] was the most common circulating genotype; however, following vaccination, G2P[4] has been more frequently detected [19,20]. In Finland, following the introduction of RotaTeq, G9P[8] and G12P[8] have now become the main genotypes, where, previously, G1P[8] dominated [21]. However, changes in genotype distribution also occurs in countries with no immunisation [22,23], so whether the vaccine directly leads to a change in genotype diversity remains unclear [24,25].

Surveillance of rotavirus genotypes has been recommended by the WHO in countries with immunisation programmes to detect and monitor strain variation and ensure vaccine effectiveness is maintained [26]. The surveillance network, EuroRotaNet, has been monitoring rotavirus diversity in 12 European countries and has reported an increase in diversity since vaccination [17,27]. As Ireland is not currently part of any European or global surveillance network, we aim to fill that current gap of knowledge.

The purpose of this study is two-fold; firstly, to report on the viruses causing gastroenteritis, including rotavirus, 2 years prior and 3 years post implementation of the Rotarix immunisation programme, and, secondly, to describe the diversity of rotavirus genotypes in Ireland.

#### **2. Results**

#### *2.1. Sample Demographics*

Ireland has a population of 4.8 million with 36% of people living in the eastern health region, which includes Dublin, the surrounding areas, and the country's largest children's hospitals [28]. The National Virus Reference Laboratory (NVRL), Dublin, provides a diagnostic and reference service for all health care regions, though testing is also provided in regional hospitals.

This study analyzed the results from pediatric (≤5 years) patient samples received at the NVRL between 1 January 2015 and 31 December 2019 for the investigation of viral gastroenteritis. In total, 11,800 faecal samples were included in the analysis, 5267 (45%) from females, 6511 (55%) from males, and 22 (0.2%) for which details were not provided. Samples tested were predominantly from the eastern health region 10,644/11,800 (90%), and of these 10,180/10,644 (96%) were from a children's hospital. Other samples were from the northern 672/11,800 (6%), western 204/11,800 (2%), midlands 139/11,800

(1%), and southern health regions 141/11,800 (1%). As vaccine history was not available for each patient, cohorts are described as vaccine-eligible, using age as a proxy for vaccination status.

During 2015 to 2019, there were 312,013 births recorded in Ireland; 159,821 males (51.2%) and 152,192 females (48.8%). To establish how representative the samples tested were, the percentage of the annual birth cohort investigated for the detection of viruses causing gastroenteritis was calculated for those aged 0–12 months in each year. In 2015, 2280/65,536 (3.5%) were tested, in 2016 2065/63,841 (3.3%) were tested, in 2017 760/61,824 (1.2%) were tested, in 2018 587/61,016 (1%) were tested, and in 2019 608/59,796 (1%) were tested.

#### *2.2. Detection of Viral Pathogens*

The most frequently detected viral pathogen in 2015 and 2016 was rotavirus, followed by norovirus. Norovirus has been detected in approximately 12% of samples each year, whereas enteric adenovirus (adenovirus subgenus F), sapovirus, and astrovirus were detected in 2.8–6.3% of samples from 2015–2019 (Table 1). The number of samples with no virus detected ranged from 51.4% to 65.0%, depending on the year.



<sup>a</sup> 185 dual infections, 10 triple infections <sup>b</sup> 154 dual infections, 11 triple infections, 1 quadruple infection, <sup>c</sup> 93 dual infections, 8 triple infections, 2 quadruple infections <sup>d</sup> 41 dual infections <sup>e</sup> 25 dual infections. Additional viruses detected in Rotarix samples, 2017: norovirus n = 3, adenovirus F n = 1, astrovirus n = 2; 2018: norovirus n = 8, adenovirus F n = 1; 2019: norovirus n = 3, sapovirus n = 1, astrovirus n = 1.

There were 1753 samples tested from vaccine-eligible children in 2017–2019, and of these 43 (2.5%) had wild-type rotavirus, 179 (10.2%) had Rotarix, 257 (14.7%) had norovirus, 113 (6.4%) had adenovirus F, 97 (5.5%) had sapovirus, and 95 (5.4%) had astrovirus detected. In this group, there were 70 dual infections and 1039 (59.3%) samples had no detectable virus.

#### *2.3. Detection of Wild-Type Rotavirus*

The median age of those testing positive for rotavirus in the pre-vaccine era, 2015–2016 (n = 1181), was significantly lower at 1.19 years (interquartile range (IQR) 0.64–1.85), compared to the median in the entire 3 years post-vaccine, 2017–2019 (n = 373) at 1.85 years (IQR 1.12–2.84) *p* < 0.0001 (Table 2).

In the 2015/16 pre-vaccine era, a total of 485/4345 (11.2%) children aged 0–1 year had detectable wild-type rotavirus. This compares with 12/1195 (1.0%) in the post-vaccine 2018 and 2019 era, representing a 91.1% relative decrease in the number of wild-type rotavirus detected in this age range. The 1–2-year age group showed a relative reduction of 79.1% when 2015/16 was compared with 2018/19; 444/1691 (26.3%) compared to 28/505 (5.5%), respectively. This contrasts with the 5–6-year age group, which showed an increase in the detection of rotavirus from 20/324 (6.2%) to 12/96 (12.5%).


**Table 2.** Number of wild-type rotavirus positive cases by age group in the pre-vaccine years (2015–2016) compared to the post-vaccine years (2017–2019).

Interquartile range (IQR).

#### *2.4. Seasonal Variation of Wild-Type Rotavirus*

Prior to vaccination, rotavirus was a seasonal infection. In 2015, the season ran from weeks 1–29, peaking in week 11; in 2016 from weeks 11–27, peaking in week 19; and in 2017 from weeks 2–30, peaking in week 11. However, in 2018 and 2019, there was no clear seasonal onset and end, and rotavirus was most frequently detected in weeks 14 and 22, respectively (Figure 1).

**Figure 1.** Number of wild-type rotavirus (Wt-RV) cases detected by month, 2015–2019.

#### *2.5. Detection of Vaccine-Derived Rotavirus (Rotarix)*

Of all 3814 samples tested in 2017–2019, 1753 (46.0%) were eligible for the vaccine, and Rotarix was detected in 179/1753 (10.2%) of these (Table 1). In addition, one sample from a vaccine eligible patient was received and tested in December 2016 and found to be positive for Rotarix. In 20/180 (11.1%) of Rotarix-positive samples, another virus was detected, most commonly norovirus.

The age at which Rotarix was most frequently detected was 2 months (Table 3). Rotarix was not detected in any samples from patients older than 7 months of age.



#### *2.6. Distribution of Genotypes in Ireland*

In total, 786/1554 (51%) samples with detectable wild-type rotavirus were genotyped. Of these, 728 (93%) were from the eastern health board, 33 (4%) northern, 14 (2%) western, 8 (1%) southern, and 3 (0.4%) from the midlands. No significant correlation was observed between genotype and region or age (data not shown).

As the total numbers of rotavirus cases decreased following the introduction of immunization in December 2016, the proportion of samples genotyped was increased to reliably detect significant changes. In 2015, 293/662 (44%), and in 2016, 242/519 (47%) positive samples were genotyped, while in 2017, 135/250 (54%), in 2018, 48/53 (91%), and in 2019, 68/70 (97%) were genotyped.

#### *2.7. Comparison of the Genotype Diversity Pre- and Post-Vaccine*

G1P[8] was the most common genotype detected in 2015, 2016, and 2017 (Table 4). Conversely, G2P[4] was the most frequently detected genotype in 2019 (Figure 2). G3P[4], G8P[8], G9P[4], G12P[6], and G2P[8] remain uncommon genotypes in Ireland, detected in five, two, four, two, and one samples, respectively, over the 5-year period.

**Table 4.** Comparison of genotype diversity between the pre-vaccine (2015–2016) and post-vaccine (2017–2019) eras. Confidence interval (CI) significance 0.95. Comparison of the genotype proportion between pre- and post-vaccine year groups by Chi-square *p* < 0.05.


**Figure 2.** Genotype diversity in Ireland 2015–2019. Data are presented as the proportion (%) of a specific genotype compared to the total genotype results. Uncommon genotypes: <1% of total results, 2015: G9P[4] n = 2; 2016: G12P[6] n = 2; 2017: G8P[8] n = 1. 2018; G3P[4] n = 1; 2019 G2P[8] n = 1, G3P[4] n = 4, G8P[8] n = 1, G9P[4] n = 2. Mixed genotypes: those with >1 G or P-type, 2015: G1/4P[8] n = 7; 2016: G8/12P[8] n = 1, G2/3P[8] n = 1, G2/9 P[4/8] n = 1; 2017: G1/3P[8] n = 1; 2018: G1/3P[8] n = 1; 2019: G8/12P[8] n = 2, G3/12P[8] n = 1, G8P[8] n = 1, G9/12P[8] n = 1. Untypable results are those where either the G or P type was untypable.

#### *2.8. Detection of Human and Equine-Like Rotavirus G3*

G3P[8] was detected in 19/535 (4%) of genotyped samples in 2015–2016, compared to 40/251 (16%) in 2017–2019, a significant increase (p < 0.0001), whilst the uncommon G3P[4] was only detected in the post-vaccine era. Four G3 strains were detected as a mixed infection. Of the 68 G3 types (63 P[8] and 5 P[4]), 17 were selected for sequencing of the VP7 gene, which identified two G3P[8] samples from 2018, containing viruses that clustered within the equine-like G3 lineage and the remaining 15 G3 samples clustered within the human lineage (Figure 3).

**Figure 3.** Phylogenetic tree of VP7 G3 rotavirus gene sequences. The tree was constructed by using the maximum likelihood method and the Tamura-Nei model [29]. Bootstrap values (1000 replicates) above 75% are shown. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. This analysis shows 46 nucleotide sequences; 17 from Irish strains identified in this study (colour coded with green circle) and 29 from reference strains in GenBank. Phylogenetic analyses were conducted in MEGA X [30].

#### *2.9. Genotypes Detected in Rotavirus Positive Samples from Those of Vaccine-Eligible Age*

There are six common genotypes circulating in Europe, namely, G1P[8], G2P[4], G3P[8], G4P[8], G9P[8], and G12P[8] [17], and all other genotypes are considered uncommon. Of the 43 samples with detectable wild-type rotavirus and of vaccine-eligible age, 37 were genotyped by RT-PCR (Table 5). Of these, 30/37 (81.1%) had a common genotype, 4/37 (10.8%) had an uncommon genotype, 2/37 (5.4%) had a mixed infection, and one sample (2.7%) could not be fully genotyped. Ten of the 37 samples (27.0%) had G3 genotype detected (seven P[8] and three P[4]).


**Table 5.** Wild-type rotavirus genotypes detected in the age group eligible for the vaccine. Six additional samples had detectable wild-type rotavirus but were unavailable for genotyping.

#### **3. Discussion**

This study describes the reduction in rotavirus detection following implementation of a national immunization program for all children in Ireland, as well as previously unknown data regarding the extent of genotype diversity during 2015–2019. Our study shows that the largest reduction in the detection of rotavirus occurred in those aged 0–12 months, where a relative decrease of 91% was achieved between 2015/16 and 2018/19. Although the vaccine status was unknown in detail, the effectiveness of the vaccination program has been clearly shown. Our results support the national data collated by the Health Protection Surveillance Centre (HPSC), where a crude incidence rate (CIR) of rotavirus for all age groups was 13.3 per 100,000 population in 2018, representing a decrease of 76%, compared to the mean CIR during 2008–2017 of 55.5 per 100,000 [5]. In addition, since the introduction of the vaccine, there was a reduction in visits to three large pediatric emergency departments with acute gastroenteritis, where median weekly presentations in 2017–2018 (126; interquartile range (IQR), 103–165) were lower than in 2012–2016 (160; IQR 128–214) (*p* < 0.001) [31]. Furthermore, an 86% (95% CI 79.3–90.2%) decrease in hospitalizations due to rotavirus has been reported nationally in those aged <1 year [32]. In our study, we found that the median age of wild-type rotavirus infection significantly increased in the years following vaccination, from 1.2 years in 2015 to 2.9 years in 2019 (*p* < 0.0001). This is consistent with the findings of other researchers, who also noted the later age of infection in the post-vaccine era [11,33]. In Ireland, vaccination uptake is recorded by individual General Practitioners and health care professionals which are submitted to the HPSC on a quarterly basis. Vaccine uptake data from Q1 2017 to Q3 2017 was unavailable, but the evidence suggests that this must have been suboptimal as there was little change in rotavirus detection in those aged 0–1 year in 2016 compared to 2017 (9.7% versus 7.5%, respectively). The substantial increase in rotavirus infection in the post-vaccine years in the 5–6-year age group who would not be eligible for the vaccine was also somewhat surprising. Although the number of children tested in this age group were lower in the post-vaccine compared to pre-vaccine years, the proportion of positives was almost double (6.2% versus 12.5%). The short timeframe is a limitation of this study; however, collection of data is ongoing, and it will be of interest to follow up on the impact of vaccination on rotavirus detection in all age groups in the post-vaccine era. A further finding in our data is a diminution of the characteristic rotavirus seasonal pattern, a phenomenon that has been noted by others following introduction of the rotavirus vaccine [17,34].

The live-attenuated vaccine Rotarix replicates in the gut of the recipient and is excreted, albeit at lower amounts compared to a wild-type infection [35]. We detected Rotarix in 10% of patients who were of vaccine-eligible age and, as rotavirus is notifiable in Ireland, this highlights the importance of differentiating between wild-type and vaccine-derived viruses, particularly when screening with a sensitive method, such as RT-PCR. By not excluding vaccine-derived rotavirus from diagnostic tests, there may be an over-estimation of rotavirus disease burden and unnecessary clinical intervention [36–38]. We identified 180 samples with detectable Rotarix, 20 (11%) of which had another virus detected, the most common being norovirus. We found norovirus to be the second most common pathogen detected after rotavirus in 2015/16, which then became the most common cause of viral gastroenteritis in our study group in the post-vaccine era. Our results are consistent with that observed in earlier studies, where norovirus is now the leading cause of viral gastroenteritis in those vaccinated for rotavirus [39–41]. Of note, sapovirus, astrovirus, and enteric adenovirus were detected in similar proportions over the 5-year time period and demonstrated no increase or decrease in detection rates following Rotarix introduction. Depending on the year, we report that 51–65% had no detectable viral pathogen. This apparent diagnostic gap highlights a further limitation of this study, in that it is quite possible that parallel samples were sent for the investigation of bacterial or parasitic pathogens, which are common causes of gastroenteritis [42,43]. Unfortunately, we did not have access to these results. In addition, other viruses, such as bocavirus, enterovirus, and parechoviruses, which may cause gastroenteritis, would not have been detected by our routine screening test.

Prior to the introduction of Rotarix, we found the circulating genotypes in Ireland were comparable to other European countries, with G1P[8] being the most commonly detected. The findings of the current study are consistent with those observed in several earlier reports from samples tested in Ireland from 1995 to 2009, where it was reported that the most commonly detected genotype was G1P[8], with fluctuating levels of G2P[4], G3P[8], G4P[8], and G9P[8] [44–49]. The current study matches those findings. However, we can report that the diversity of genotypes increased in the years following the introduction of a vaccine and that, in 2018/19, G1P[8] was no longer the most common genotype. Furthermore, genotypes detected in children eligible for the vaccine was more varied than those detected in the vaccine-ineligible cohort. With regards to wild-type and vaccine rotavirus strains, they can be described in terms of being homotypic, partly heterotypic, and fully heterotypic based on the G and P proteins. For instance, the monovalent G1P[8] Rotarix vaccine is homotypic to other circulating G1P[8] strains (both proteins are the same), G12P[8] is partly heterotypic (one protein different), and G2P[4] is fully heterotypic (both G and P proteins are different) [50]. Rotarix provides exposure to G1P[8] rotavirus among infants, with protection that is likely to be higher against homotypic strains than heterotypic strains, such as G2P[4]. This suggests that natural infection leading to disease is more likely to be caused by such heterotypic strains [19,51] and that a vaccinated population could possibly drive selective pressure, increasing the likelihood of these genotypes to circulate in the community [52,53]. That being said, the monovalent vaccine Rotarix provides significant protection from G1, G2, G3, G4, and G9, and efficacy against severe G2 rotavirus gastroenteritis was as high as for other rotavirus types [54]. Clearly, the immune response to rotavirus infection is a complex issue, with a previous report suggesting that type-specific neutralizing antibodies induced by the vaccine against VP7/VP4 epitopes are not solely responsible for a protective effect [55]. The report proposes that, as there are a limited number of diverse circulating strains worldwide, these antibodies are not driving long-term selective pressure, which itself would favor antigenic drift or the emergence of novel genotypes.

Interestingly, the genotype G12P[8] was not circulating widely at the time of Rotarix and RotaTeq vaccine development; however, it has now become established as an increasingly common genotype [17]. A large study in the USA found G12P[8] more frequently than any other rotavirus genotype in fully vaccinated children [56]. Another example of an uncommon genotype becoming more prevalent is the recently emerged equine-like G3 strain, first identified in Japan in 2013 [57] but now detected world-wide [58–60]. We identified, for the first time in Ireland, two samples from 2018 that clustered within the equine-like G3 lineage. A further 15 G3 samples were sequenced and all clustered in the human lineage, suggesting that the equine-like lineage has not yet become established in Ireland compared to other countries [59,61]. Of note, five of the uncommon human G3P[4] strains were

also identified in our study in 2018–2019, and this strain has been detected before in Ireland in 2006/07 [48]. The detection of uncommon genotypes, along with the additional potential for zoonotic reassortment [62,63], reinforces the WHO recommendation for surveillance, emphasising the need for continued monitoring of rotavirus vaccine efficacy against emerging rotavirus.

Several important limitations need to be considered for our study group. Firstly, the results are somewhat biased due to the observational nature of the study and samples tested would have been from those with moderate to severe gastroenteritis that warranted clinical investigation. In addition, there was no denominator for the population not suffering from symptoms of viral gastroenteritis, so we are unable to calculate incidence and prevalence of rotavirus infection. Furthermore, with no access to vaccination data, it was not possible to determine vaccine effectiveness rates or describe definitive vaccine failures. However, due to the large data set (n = 11,800), we can show relative reductions in the detection of rotavirus and the changes in the diversity of circulating genotypes. The geographical distribution of samples is nationwide, although the data are skewed to some extent due to the density of the population in Dublin and the location of the main children's hospitals, and therefore samples were predominantly from the eastern health region. It should also be noted that the number of samples tested has decreased year on year. The reason for this may be the decrease in symptomatic children or the increase in localized testing and possible availability of point of care assays. Finally, it was not possible to categorize samples as the community or hospital acquired and we could not identify samples belonging to outbreaks.

We have shown that rotavirus continues to circulate in the pediatric population, albeit in low numbers, and this is expected to decrease further with the increasing cohort of vaccinated children. Binary genotypic classification is useful to establish circulating genotypes and can be used for reassortment studies of the VP7 and VP4 encoding genes; however, whole genome genotyping is required for a more detailed analysis of the virus. Indeed, a future aim from this ongoing study is to perform whole genome sequencing from samples in this dataset to allow identification of possible reassortment of non VP7/VP4 genes or mutation events. In addition, all samples identified as G3 will be categorized as either equine-like or of human lineage. By collaborating with clinicians at the children's hospitals, it is hoped that any sample from a child with rotavirus with a full or partial vaccine history will be referred to the NVRL for whole genome sequencing to establish definitive strains circulating in this group of children.

In conclusion, we describe the detection and characterization of rotavirus in pediatric samples circulating in Ireland over a 5-year time period. We show that, following the introduction of Rotarix, there is a relative reduction in the number of rotavirus infections diagnosed, coinciding with an increase in genotype diversity, along with the first recorded detection of an equine-like G3 strain in Ireland.

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

#### *4.1. Study Design*

This opportunistic study presents the results of faecal samples from pediatric samples (≤5 years) investigated for viral gastroenteritis at the National Virus Reference Laboratory (NVRL), Dublin, Ireland. Test results for wild-type rotavirus, vaccine-derived rotavirus (Rotarix), norovirus, sapovirus, astrovirus, and enteric adenovirus subgenus F were obtained with genotype and sequence results, if available. Samples dated 1 January 2015 to 31 December 2019 were included in the study. Samples dated 1 January 2015 to 31 December 2016 were designated as "pre-vaccine". The first doses of Rotarix were given for those aged 2 months from the 1 December 2016. Only one sample was received from a 2-month old in December 2016, and this patient had detectable Rotarix. Samples dated 1 January 2017 to 31 December 2019 were designated as "post-vaccine". Routine testing for Rotarix was introduced into the NVRL from 11 December 2017. Rotavirus-positive samples received from 1 December 2016 to 10 December 2017 were tested for Rotarix, retrospectively.

#### *4.2. Annual Birth Cohort in Ireland*

The Central Office of Statistics provides the annual number of births in Ireland [64]. The number of annual births by year are: 2015: 65,536 (33,480 males, 32,056 females); 2016: 63,841 (32,709 males, 31,132 females); 2017: 61,824 (31,779 males, 30,045 females); 2018: 61,016 (31,298 males, 29,718 females); 2019: 59,796 (30,555 males, 29,241 females,). Data for 2018 and 2019 are provisional. The overall male: female ratio for 2015–2019 was 1.1: 1 (51.2% versus 48.8%).

#### *4.3. Data Analysis*

Data were extracted from the NVRL Laboratory Information Management System and analyzed in Excel. All samples were assumed to be from symptomatic patients. Samples with no date of birth recorded or duplicate samples were excluded from the study. Patients were de-identified, and the variables recorded in the database were patients' age at sample collection, sex, sample date, geographical region, and test result(s). Geographical regions were categorized as eastern (which includes Dublin), western, southern (south and south-east), northern (north-west and north-east), and midlands (midlands and mid-west), as defined by the Health Service Executive areas used by the Health Protection Surveillance Centre [5].

#### *4.4. Sampling Strategy for Genotyping and Sequencing of Samples*

To determine the sample size required to reliably detect a change in genotype frequency, a sample size calculator was used (CL95%, www.openepi.com), and then a random selection of wild-type rotavirus samples was selected for genotyping. In addition, all uncommon genotypes and a random subset of genotyped samples was selected for sequencing of VP7 and VP4 genes. A subset of those identified as G3 and were sequenced and analyzed to determine a human or equine-like lineage.

#### *4.5. Seasonality*

Seasonal onset, peak, and end were calculated: Onset: First of 2 consecutive weeks, where the median percentage of positive results was >10%. Peak: Week with the highest proportion of positive samples. End: Last of two consecutive weeks, where the median percentage was <10%. Denominator: total samples tested; numerator: the number of positive samples.

#### *4.6. Vaccine Eligibility*

The vaccine status was unknown, and patients were categorised by vaccine eligibility. Vaccine-eligible samples were those born after 1 October 2016 and were ≥2 months of age. Vaccine-ineligible samples were those born after 1 October 2016 and were <2 months of age or were born prior to 1 October 2016 and were aged 0–5 years of age.

#### *4.7. Laboratory Methods*

Upon receipt into the laboratory, approximately 20% *w*/*v* suspension of the fecal sample was prepared in 400 µL Stool Transport and Recovery Buffer (Roche) and 400 µL external lysis buffer (Roche). A total of 450 µL of the suspension was extracted by Roche MagNAPure 96 and eluted into 100 µL. During extraction, Brome Mosaic Virus RNA (University of Indiana) was added as an internal control (IC) at 1 pg/µL to the sample prior to extraction. The eluates were tested in five one-step RT-PCR assays, as previously described [65–70]. Briefly, eluates were tested in a 25 µL or 10 µL reaction mixture (depending on the 96- or 384-well format, respectively), containing 2× Superscript™ III Platinum One-Step qRT-PCR mix (Invitrogen), as per product insert. Final concentrations of primers and probes ranged from 80 nM to 400 nM, depending on the target. Each sample eluate was tested by five (RT-)PCR reactions, namely: (i) norovirus G1/G2/IC; (ii) adenovirus F/pan-rotavirus/IC; (iii) Rotarix; iv) astrovirus/IC; v) sapovirus/IC. Amplification was performed on the ABI 7500 Fast (96-well format) or the ABI Viia7 (384-well format) instrument under the following conditions: 15 mins 50 ◦C, 2 mins 95 ◦C, 38 cycles of 15 secs 95 and 30 secs 60 ◦C (56 ◦C for norovirus). Amplification data was collected and analyzed with Sequence Detection Software version 2.3 or the Viia7 software version 1.2.1 (both from Applied Biosystems).

Genotyping was either by a multiplex RT-PCR [71,72] or by hemi-nested RT-PCR, as described previously [73], with fragment visualization and size determination performed on the TapeStation (Agilent software, version 2200). Samples with an indeterminate G or P type were tested by both methods before being categorized as untypable. A selection of genotypes were confirmed by Sanger sequencing of the VP7 and VP4 genes, using previously described methods [73] on the ABI 3500Dx genetic analyzer (Applied Biosystems) and typed using the RotaC typing tool [74] or by the Basic Local Alignment and Search Tool, BLAST (http://www.ncbi.nlm.nih.gov/blast/Blast.cgi). G3 VP7 sequences (450 nucleotides) were aligned with appropriate reference sequences using ClustalW. Phylogenetic analyses were conducted in MEGA X [30] using the maximum likelihood method, with 1000 bootstrap replicates, based on the Tamura-Nei model [29]. This model was selected as it generated the lowest Bayesian information criterion (BIC) score in MEGA X.

#### *4.8. GenBank Accession Numbers*

Partial VP7 fragments of the equine-like G3 strains identified in this study were deposited in GenBank under the following accession numbers: strains; MT475885 and MT4758866, whereas the human-lineage G3 VP7 fragments were MT537569-537583.

#### *4.9. Statistical Analysis*

The study was observational and therefore most data presented was descriptive. The median age of rotavirus infection in the pre- and post-vaccination groups was compared using Mann–Whitney *U* test. The Chi-square test for proportions was used to compare genotypes in the pre- and post-vaccination groups. *P* values for both tests of ≤0.05 were considered statistically significant. Confidence intervals (95%) were calculated using the Wilson method for a proportion of the genotypes detected. Statistical analysis was performed using SPSS 26 (IBM Corp; Armonk, NY, USA) software or www.openepi.com.

#### *4.10. Ethical Statement*

All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or research committee and with the 1964 Helsinki declaration and its later amendments or comparable standards. This study was approved for ethical exemption by University College, Dublin LS-E-17-09.

**Author Contributions:** Conceptualization, Z.Y., S.C., and C.F.D.G.; data curation, S.C.; formal analysis, Z.Y.; investigation, Z.Y. and S.C.; methodology, Z.Y., S.C., and C.F.D.G.; project administration, Z.Y. and S.C.; resources, Z.Y., G.T., A.C., and C.F.D.G.; supervision, S.C. and C.F.D.G.; validation, Z.Y., S.C., and C.F.D.G.; visualization, Z.Y., S.C., and C.F.D.G.; writing—original draft, Z.Y.; Writing—review and editing, Z.Y., S.C., J.D., G.T., and C.F.D.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding

**Acknowledgments:** The authors would like to thank all clinicians and general practitioners for referring samples to the National Virus Reference Laboratory, Dublin, and to acknowledge the hard work of all the laboratory staff for processing and testing samples for the gastroenteritis screen.

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

#### **References**

1. Hungerford, D.; Smith, K.; Tucker, A.; Iturriza-Gomara, M.; Vivancos, R.; McLeonard, C.; N, A.C.; French, N. Population effectiveness of the pentavalent and monovalent rotavirus vaccines: A systematic review and meta-analysis of observational studies. *BMC Infect. Dis.* **2017**, *17*, 569. [CrossRef] [PubMed]


© 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/).
