Impact of Vaccination on Rotavirus Genotype Diversity: A Nearly Two-Decade-Long Epidemiological Study before and after Rotavirus Vaccine Introduction in Sicily, Italy

Sicily was the first Italian region to introduce rotavirus (RV) vaccination with the monovalent G1P[8] vaccine Rotarix® in May 2012. In this study, the seasonal distribution and molecular characterization of RV strains detected over 19 years were compared to understand the effect of Rotarix® on the evolutionary dynamics of human RVs. A total of 7846 stool samples collected from children < 5 years of age, hospitalized with acute gastroenteritis, were tested for RV detection and genotyping. Since 2013, vaccine coverage has progressively increased, while the RV prevalence decreased from 36.1% to 13.3% with a loss of seasonality. The local distribution of RV genotypes changed over the time possibly due to vaccine introduction, with a drastic reduction in G1P[8] strains replaced by common and novel emerging RV strains, such as equine-like G3P[8] in the 2018–2019 season. Comparison of VP7 and VP4 amino acid (aa) sequences with the cognate genes of Rotarix® and RotaTeq® vaccine strains showed specific aa changes in the antigenic epitopes of VP7 and of the VP8* portion of VP4 of the Italian RV strains. Molecular epidemiological surveillance data are required to monitor the emergence of novel RV strains and ascertain if these strains may affect the efficacy of RV vaccines.


Introduction
Before the introduction of rotavirus (RV) vaccines, RV acute gastroenteritis (AGE) in infants and children under five years of age led to approximately 500,000 deaths each year, mainly in developing countries, where the access to safe drinking water and sanitary care are sub-optimal. With the introduction of RV vaccines, the current annual mortality estimates in children under five years of age range from 122,322 to 21,575, with deaths still mainly occurring in developing countries [1]. In developed countries before RV vaccine introduction, RV infections resulted in a substantial economic impact on healthcare systems and families due to medical visits, emergency room access and hospitalizations [2,3]. In 2006 and 2008, two oral RV vaccines, a pentavalent human-bovine reassortant vaccine (RotaTeq ® [Merck & Co] West Point, PA, USA, RV5) and a monovalent vaccine (Rotarix ® [GSK], Rixensart, Belgium, RV1), were licensed [4]. Both vaccines demonstrated high efficacy (>85%) against severe AGE episodes [5,6]. To date, almost 100 countries have introduced RV vaccines, including countries in Sub-Saharan Africa, the Americas, Europe, and the Eastern Mediterranean/Middle East region (The RV Organization of Technical Allies (ROTA) Council. http://rotacouncil.org/vaccine-introduction/global-introductionstatus/ accessed on 7 July 2021). The impact on severe RV and all-cause diarrhea has been notable in countries that have introduced the vaccine. Epidemiological surveys on RV prevalence and molecular characterization of the viral strains circulating in the pediatric population are important to monitor vaccine effectiveness and to promptly detect novel RV strains/variants that might emerge under vaccine pressure and potentially affect vaccine-acquired protection.
Sicily was the first Italian region to introduce RV vaccination into the routine childhood immunization schedule (June 2012). Since then, Rotarix ® (monovalent G1P [8] (Glaxo-SmithKline [GSK]Vaccines, Rixensart, Belgium) has been the vaccine offered by the local healthcare authority, being administered on a two-dose schedule at 3 and 7 months of age. Vaccine coverage in 2012 was very low (6%) in Sicily but progressively increased from 20% in 2013 to 54% in 2019. In countries implementing universal RV vaccination, surveillance data have shown that vaccine introduction caused a reduction in RV disease burden, but also an increase in selective pressure exerted on strains circulating into human populations [14][15][16].
In Palermo, Sicily, more than 30 years (uninterruptedly since 1985) of RV surveillance data are available. This represented a unique opportunity to monitor changes in the prevalence of RV genotypes and to assess temporal patterns of variation over a large time window and a broad number of samples. This remarkable epidemiological setting provides a robust context through which to evaluate the local effects of the introduction of RV vaccination.
In this study, we compared the incidence of severe RV AGE, the seasonal distribution of RV genotypes, and the molecular epidemiological data of RVs causing infections in Sicily before (January 2002 to December 2012) and after (January 2013 to December 2020) the introduction of the RV vaccine.

Seasonality and G/P Temporal Distribution
In total, 7846 samples collected during 19 consecutive years of hospital-based surveillance, from January 2002 to December 2020, were retrospectively analyzed to evaluate the epidemiological changes in RV ecology after the introduction of the monovalent G1P [8] vaccine (Rotarix ® , RV1) in Sicily, Italy. The overall rate of RV-positive samples detected during the entire study period was 22.6%. Out of the 1773 RV-positive samples detected over the whole study period, 1157 (65.2%) were detected from 2002 to 2012, and 616 (34.7%) from 2013 to 2020. All the RV-positive samples were used for epidemiological analyses of the pre-and post-vaccine era, respectively (Table 1). RV infection yearly prevalence rates ranged from 15.54% to 54.21% in the pre-vaccine period, with a combined pre-vaccine prevalence of 36.18%, and from 6.20% to 26.59% after vaccine introduction, with a combined prevalence of 13.36% over the post-vaccine period. A significant difference in pre-and post-vaccine periods was observed when comparing the yearly prevalence of RV infection in the two study periods (two-sample t-test; p < 0.05). In pre-vaccine era, from 2002 to 2012, RV infection exhibited increased seasonal activity from the end of winter to the end of spring, with prevalence ranging from 56.96% to 44.24% from February to May. While in the post-vaccine era, from 2013 to 2020, no period of specially increased circulation could be observed; only a plateau of prevalence rates above or approaching 20% (range: 16.98-21.83%) extending from the end of spring to summer was observed ( Figure 1). Interestingly, when comparing data from pre-and post-vaccine periods, a 1.2-to 6.5-fold decline in monthly RV infection rates was observed after vaccine introduction ( Table 2).
Among the 1773 RV-positive samples detected during the study period, an RV genotype was determined for 1651 (93.1%) samples. Specifically, 1082 out of 1157 (93.5%) RVpositive specimens were genotyped in the pre-vaccine period and 569 out of 616 (92.37%) in the post-vaccine period.

VP7 and VP4 Phylogenetic Analyses
For a selection of RVs belonging to different G-and P-types, the VP7 and VP4 genes were sequenced and analyzed phylogenetically (Supplementary Figure S1). Among the G1P [8] strains, a single sub-lineage, Ic, which had been circulating since 2004, was predominant (75% of the sequenced strains) until 2019, when it was replaced by sub-lineage IIa (Figure 3a). In the VP7 gene, the Italian G1-Ic RVs showed a percentage of nucleotide (nt) and amino acid (aa) identity ranging from 95 to 100% similarity to each other and from 91 to 95% nt and 91 to 97% aa similarity to the Rotarix ® strain (sub-lineage IIa), and ranging from 89 to 91% nt and from 90 to 94% aa similarity to the RotaTeq ® strain (lineage III). All the G1-IIa strains detected were vaccine strains showing a 98.3-100% nt and 98-100% aa identity similarity to the Rotarix ® strain. The Italian G2P [4] RVs detected from 2003 to 2019 clustered into lineage IV, in two distinct sub-lineages, a-1 and a-3, alternating over time ( Figure 3b). The two sub-lineages showed an intra sub-lineage, nt, and aa identity of 96 to 100% and an inter sub-lineage nt and aa identity of 91-97% and 96-99%, respectively. The VP7 of the RotaTeq ® G2 strain clustered into lineage II and showed a 92-94% nt and 92-97% aa identity similarity to the Italian G2-a-1 strains and a 91-94% of nt and 91-99% aa identity similarity to GII-a-3 strains. Among the G3P [8] RVs, the Italian strains detected from 2004 to 2018 clustered in lineage III together with "classical" human strains (96.5-100% nt identity and 95.8-100% aa identity), while the G3P [8] strains that had been detected since 2019 segregated into lineage I together with novel equine-like RVs (98.8-100% nt identity and 97.1-100% aa identity). A contemporaneous circulation of G3P [8] equine-like strains, clustering into lineage I (98-100% nt and a 96-100% aa identity), and "classical" G3P [8] strains of lineage III (97-100% nt and a 98-100% aa identity) was observed in 2020 ( Figure 3c).
Sequence analyses of the VP7 gene of G4P [8] RVs showed that all Italian strains detected from 2002 to 2020 fell together into lineage Ic, with a 98-100% nt and a 96-100% aa similarity to each other. These G4 strains showed a 94-96% nt and 86-90% aa identity similarity to the RotaTeq ® G4 strain. G9P [8] RVs were frequently detected in both the preand post-vaccine periods, and they all belonged to a major sub-cluster within lineage III (94-100% nt and 90-100% aa identity). All the Italian G12P [8] RVs clustered into lineage III, but into two different sub-lineages (III-a and III-b). In particular, the G12 strains detected in 2012 and 2014 were equally distributed into sub-lineages III-a and III-b, whilst all 2015 and 2016 G12 strains clustered in sub-lineage III-b, and the 2020 G12 strains were segregated into sub-lineage III-a. Within these sub-lineages, the identity was 96-100% nt and 98-100% aa for lineage III-a and 95-100% nt and 95-100% aa for lineage III-b (Figure 3d).  In the VP8* portion of VP4, all the VP4 P [8] strains, regardless of the associated G-type (G1, G3, G4, G9, and G12), clustered into lineage III, showing 92 to 100% similarity in nt and aa identities to each other. When compared to the vaccine strains, the P [8] strains showed 87 to 91% nt identity and 86 to 92% aa identity similarity to the Rotarix ® strain (Lineage I), and 89-94% nt and 88-96% aa identity similarity to RotaTeq ® strain (lineage II). The Italian P [4] strains clustered into lineage IV (93-100% of nt and 91-100% aa identities). Phylogenetic analysis showed a lower genetic diversity without any apparent temporal pattern of distribution among P[4] RVs.
The letters represent the peculiar amino acid substitutions of neutralization epitopes in the VP8* portion of VP4 compared with homologous sequences of the Rotarix ® and RotaTeq ® vaccine strains, were showed in Table 4.

Discussion
Then RV vaccine was introduced in Sicily, Italy, in June 2012. The present study investigated the epidemiology of RV infection before (from 2002 to 2012) and after (from 2013 to 2020) vaccine introduction. Based on the epidemiological data, the increasing vaccine coverage induced a marked decrease in RV infection prevalence from 2013 onward, with an overall prevalence rate declining from 36.18% in the pre-vaccine period to 13.36% in the post-vaccine period. Moreover, we observed the loss of the previously strong winter/spring seasonality of RV infection, with a remarkable decrease (ranging from 1.2 to 6.5-fold) in RV infections during winter and spring in the post-vaccine period ( Figure 1, Table 2). Typically, RV infection exhibits peaks of seasonal activity during the winter season in temperate climates and throughout the year in tropical climates [17]. A delay and a blunting of seasonal peaks were also observed in other temperate climate countries where vaccination had been introduced, compared to the pre-vaccine period [18][19][20][21][22].
Based on the epidemiological information generated in our hospital-based surveillance, changes in the distribution of RV genotypes in Italy occurred after the introduction of RV vaccination. In particular, the previously predominant G1P [8] genotype, which accounted for more than 60% of RV infections in the pre-vaccine period, drastically decreased from 2014 to 2020. This strain was replaced by different genotypes in the post-vaccine period. A marked genotype diversity was observed in 2014-2015 when G1P [8] strains co-circulated at low prevalence with several other genotypes, i.e., G2P [4], G4P [8], G9P [8], and G12P [8]. The predominance of G1P [8] RV prior to vaccine introduction has been thoroughly described, as well as fluctuations in RV genotypes distribution and a switch in predominant genotypes, from G1P [8] to G2P [4] following Rotarix ® introduction [4,[23][24][25][26].
In the present study, in the post-vaccine era, the alternate prevalence of a variety of RV genotypes was observed, including G2P [4] in 2015 and 2017, G9P [8] in 2016-2017, G3P [8] in 2018-2019, and G12P [8] in 2020 (Figure 2). The emergence of these genotypes could have been driven by the force of vaccine selection, which can impose a selective pressure on circulating strains or by a selection of mutant viruses that are not effectively neutralized [4]. Alternatively, this could be the effect of natural seasonal fluctuations and global emergence of novel strains generated by natural variation or re-assortment events between human and animal strains [27,28]. In Palermo, RV surveillance activity spanned more than three decades, and the genotypes circulating in our relatively small geographic area mostly reflected the epidemiological scenario observed in Europe and in most developed countries [29,30]. Novel genotypes, such as G9P [8] in 1999-2000 and G12P [8] in 2012, emerged in Italy, acquiring epidemiological relevance over time [31,32]. Recently, novel equine-like G3P [8] strains, likely deriving from a reassortment event with equine RV strains, spread in several countries, becoming the prevalent genotype in Italy in 2018-2019 [33]. Interestingly, the circulation of equine-like G3P [8] strains abruptly dropped in 2020, replaced by G12 RV strains. Similar fluctuations of prevalent RV genotypes, with G1P [8] being replaced first by G9P [8] and more recently by G12P [8] and G3P [8] genotypes, were also observed in Brazil after the introduction of universal vaccination [34]. Several studies showed a cyclic rise in detection rates of non-G1P [8] strains in both vaccinated and non-vaccinated children [28,30,34,35]. Unfortunately, the vaccination status of the children enrolled in this study was not available and it was not possible to evaluate the differential circulation of RV types in vaccinated and non-vaccinated groups.
Genetic drift has generated several lineages and sub-lineages within RV genotypes, providing some strains with increased fitness or the ability to escape partially population immunity, and this could be related to the epidemic peaks of RV infections [15]. In previous studies, it was hypothesized that G1P [8] continuous circulation was due to the emergence end re-emergence of different lineages and sub-lineages [36]. Sub-lineage Ic, only sporadically detected in 1996-1997, 2002, and 2004, was predominant among G1P [8] strains after 2007 in Italy [29]. In this study, G1 RVs of sub-lineage Ic were found to circulate both in the pre-and post-vaccine period. A marked intra-genotype diversity, with different lineages replacing each other, was observed within G2P [4], G3P [8] and G12P [8] genotypes (Figure 3), while a single lineage of G9P [8] (lineage III) and G4P [8] (lineage Ic) RVs was detected.
The genetic diversity observed in the VP8* portion of the VP4 was lower, probably due to structural and functional constraints during attachment, penetration, and maturation of the virion. In the present study, a single VP4 lineage (P[8]-III) was observed in all P [8] strains detected over the entire study period, regardless of the associated G-type. As an exception to this rule, the G1P [8] vaccine-derived strains segregated into lineage P[8]-I.
Upon comparison of VP7 and VP4 aa sequences with the cognate genes of Rotarix ® and RotaTeq ® vaccine strains, the Italian RVs showed specific aa changes in the antigenic epitopes of VP7 (7-1a, 7-1b, and 7-2) and of the VP8* portion of VP4 (8-1 to . The VP7 sequences of Italian G1-Ic RVs, both in the pre-and post-vaccine period, showed up to four and six conserved aa substitutions with respect to Rotarix ® (A1CB052A/ 1988/G1P [8]) and RotaTeq ® (WI79-9/1992/G1P7 [5]), respectively (Table 3a). Changes at residues 94, 123, and 217 are correlated with neutralization resistance and have been also observed in G1P [8] strains circulating in several countries [37,38]. Peculiar amino acid changes, mainly located in the 7-1a and 7-b epitopes of VP7, were detected in all G2 and G4 Italian strains circulating over the entire study period. Some VP7 changes (A87T, D96N, and S213D) have been associated with the emergence of G2P [4] in different countries, following vaccine introduction, possibly due to differential antigenic pressure triggered by the monovalent G1 vaccine [15,27]. However, in our study, G2P [4] strains did not acquire particular relevance after vaccine introduction. On the other hand, a significant epidemiological change was observed in Italy in the G3P [8] genotype, with the emergence of equine-like G3 strains in 2017. Classical and equine-like G3P [8] strains showed peculiar patterns of aa substitutions with respect to the sequence of RotaTeq ® vaccine (Table 3d). Conserved aa changes, mainly located in the 8-1 and 8-3 epitopes of VP4 protein, were observed in Italian P [8] strains of lineage III, with respect to Rotarix ® (A1CB052A/1988/G1P [8]) and RotaTeq ® (WI79-4/1992/G6P [8]) vaccine strains of lineage I and II, respectively (Table 4). Such substitutions were peculiar of the P[8]-III lineage, which represents the most common VP4 type in globally circulating RVs [39].
The major antigenic changes observed in equine-like G3P [8] might allow these strains to successfully spread globally in a relatively short period [38]. The identification of notable aa substitutions in each genotype detected in both pre-and post-vaccine era might be responsible for a fitness advantage, which allowed such rotaviruses to maintain their circulation without being affected by vaccines. Otherwise, as a result of immunization strategies, the decreased circulation of RVs among humans may have reduced their innate tendency to mutate, as observed for circulating vaccine-derived polioviruses, which preferentially mutate and cause the disease in communities with low immunization rates (https://polioeradication.org/polio-today/polio-now/this-week/circulatingvaccine-derived-poliovirus/ accessed on 25 February 2022).

Study Population
To evaluate the RV variability in pre-and post-vaccine era, data collected in eight consecutive years after vaccine introduction (from January 2013 to December 2020) were compared to those collected in the 10 years before vaccine introduction (from January 2002 to December 2012). Pre-vaccine period included the entire year 2012 since vaccination was introduced only in June and coverage was considered to be too low (6%) to affect RV circulation. Overall, a total of 7846 children under 5 years of age, hospitalized with AGE at the "G. Di Cristina" Children's Hospital of Palermo, Italy, between 2002 and 2012 (pre-vaccination period, n = 3213) and from 2013 to 2020 (post-vaccination period, n = 4633), were enrolled (Table 1). This study is a retrospective study using hospital-based surveillance data and represents an extension and integration of a previous study, summarizing the results of a 27-year-long RV surveillance period in the pre-vaccination period [29].

G and P Genotyping and Phylogenetic Analyses
All stool samples were subjected to RNA extraction (QIAamp Viral RNA MiniKit, Qiagen, Hilden, Germany) and analyzed using random primers reverse transcription shortly after samples were collected [40] and real-time PCR with specific primers targeting the NSP3 gene of RV as described previously [41]. Specimens in which RV RNA was detected at a cycle threshold (Ct) value ≤ 30 were considered positive and were eligible for molecular characterization. RV-positive specimens (1495) were further analyzed to determine the G/P genotypes, according to the established binomial classification, using a mixture of specific primers (genotypes G1-G4, G6, G8, and G9-G12 and P [4], P [6], P[8]-P [11], and P [14], for VP7 and VP4 types, respectively) [42][43][44].
In addition, RV-positive fecal samples with sufficient quantity, representative of the entire study period, were submitted to sequence analyses. The nearly full-length VP7 gene and the VP8* portion of the VP4 gene were amplified with consensus primers and the sequences were determined by direct sequencing. Sequence alignment was performed with CLUSTAL W [45]. Phylogenetic analysis was carried out using the MEGA software version X [46] using the Kimura 2-parameter model as a method of substitution and the maximum likelihood method to construct phylogenetic trees from partial sequences of VP7 and VP4. Sequences were assigned to lineages/sulineages as described in previous literature [29]. Representative sequences of RVA strains belonging to different genotypes, used to generate phylogenetic trees, have been deposited in GenBank.

Statistical Analyses
The yearly pre-and post-vaccine prevalence was calculated as rotavirus-positive samples/samples tested, and a two-sample t-test was used to compare the yearly prevalence rates recorded in each of the two study periods. A p-value ≤ 0.05 was considered statistically significant. All the analyses were performed with http://vassarstats.net, accessed on 15 January 2022.

Conclusions
The introduction of anti-RV universal vaccination in Sicily encouraged monitoring of its immune-epidemiological effects. Comparative analysis of RVs detected before and after the introduction of vaccination, besides a marked reduction in the prevalence of infection, also showed increased genetic diversity in RVs. The post-vaccine emergence of G3P [8] equine-like strains and the increasing circulation of G12P [8] RVs highlight the need for molecular epidemiological surveillance in order to promptly reveal the emergence of novel RV strains, and eventually ascertain if these RV variants may escape vaccineinduced immunity.