From 2018 to 2021, a total of 356 stool samples were collected from hospitalized children with AGE. In our previous surveillance conducted in Sidoarjo, Indonesia, from 2015 to 2017, the G9P[4] strain was not detected. Interestingly, in this study, we first detected the G9P[4] strain in 2018, and G9P[4] strain continued to be detected in 2019 and 2020, albeit with low prevalence (2 out of 44, 4.5% in 2019 and 1 out of 7, 14% in 2020). Remarkably, in 2021, the majority of RVA-positive specimens (19 out of 22, 86.4%) were identified as G9P[4], while the remaining three strains (13.6%) were G9P[6].
Among them, a total of 26 samples, collected between 2018 and 2021 (Table 1) and tested positive for rotavirus kit, were analyzed by whole genome sequencing using NGS Illumina MiSeq technology. The most predominant genotype was G9P[4] (24/26, 92.3%) followed by G9P[6] (2/26, 7.7%), respectively. The genome constellation of G9P[4] (6/19, 31.5%) was G9‐P[4] (G9-P[4]-I2-R2-C2-M2-A2-N2-T2-E2-H2), i.e., the DS-1-like backbone. The reassortant strain found in G9P[4] (14/24, 58.4%), including the G9-VP7 and N1-NSP2 genes with a DS-1-like genotype constellation (G9-P[4]-I2-R2-C2-M2-A2-N1-T2-E2-H2) (Table 2). The reassortant strain found in all G9P[6] genotypes, including the G9-VP7 and N1-NSP2 genes with a DS-1-like genotype constellation (G9-P[6]-I2-R2-C2-M2-A2-N1-T2-E2-H2) (Table 2).
Phylogenetic tree analysis
To assess the evolutionary dynamics of Indonesian reassortant G9 RVA strains at the lineage level, our phylogenetic analysis incorporated sequences from the published G9P[4] RVA strains available in the DNA database. Additionally, we included a set of reference strains for genotype 2 RVA for lineage and sub-lineage classification21. Furthermore, to elucidate the relationship between the DS-1-like backbone of the newly detected reassortant G9 RVA strains in this study and the equine-like G3P[8] strains found in Indonesia during 2015-2016, we analyzed the nucleotide sequences carrying genotypes I2, R2, C2, M2, A2, N2, T2, and H2 from these equine G3P[8] DS-1-like strains (Supplementary Fig. 1).
Based on the recently proposed lineage classification for globally circulating DS-1-like strains and G9 VP7 genes13, the Indonesian reassortant G9P[4]/P[6] strains were assigned to the following significant lineages: major sub-lineage III for VP7; lineages IV for VP4, VP2, NSP1, and NSP5; lineages V for VP1, VP6, NSP2-N2, and NSP3; and lineage VI for VP3 and NSP4.
The phylogenetic trees showed that the branch topology of most segments of the Indonesian reassortant G9 RVA strains detected during 2018 to 2021 fell within the same monophyletic lineage, except for the NSP3 and NSP4 genes of three G9 strains detected in 2021 and 2022 (SOEP-733, SOEP-735 and SOEP-740). These results indicate a high nucleotide identity of 99.2% among the sequences of the Indonesian G9P[4]/P[6] strains (Fig. 1). The genetic background remains consistent among G9P[4] strains found sporadically between 2018 and 2021, when they circulated at a high prevalence.
In the VP7 phylogenetic tree, Indonesian reassortant G9 RVA strains wereclustered within the major lineage sub-lineage-III, comprising 103 G9 strains retrieved from the GenBank database spanning from 2006 to 2020 (Fig. 1A) . While the major sub-lineage III includes the G9P[4] strains identified in India, Denmark, and the USA between 2011 and 2015 as well as G9P[4] strains detected in 2018 and 2020, along with G9P[8] strains detected in Indonesia in 2018 (DSA48, DSA50, and DSA62), the newly identified G9 strains in this study exhibits a closer genetic affinity with G9P[4] strains detected in the Czech Republic in 2018 (CZE/H186/2018 and CZE/H187/2018), Russia in 2018, and Italy from 2019 to 2020. The Indonesian reassortant G9 DS-1-like strains formed a monophyletic lineage with G9P[4] strains detected between 2018 and 2020, sharing approximately 99.6% nucleotide identity within their VP7 genes (Fig. 1).
In the VP4 phylogenetic tree, all of the 23 Indonesian G9P[4] strains were clustered within lineage IV, alongside VP4 genes from previously identified G9P[4]-DS-1-like strains (Fig. 2) . Similar to the VP7 genes, the VP4 gene of Indonesian G9P[4] strains exhibited the closest relationship to the VP4 genes of G9P[4] DS-1-like strains detected in the Czech Republic in 2018 and in Italy in 2020 (ITA/PA374-20/2020, CZE/H186/2018, and CZE/H187/2018). The nucleotide similarities among them ranged from 99.6% to 99.7%. However, the VP4 gene of the reassortant G9P[6] strain detected in Indonesia in this study exhibited the highest similarity to VP4 genes from equine-like G3P[6] DS-1-like strains identified in Indonesia during 2015-2016 (Fig. 2).
In the NSP2 phylogenetic tree, 14 samples exhibited reassortment in N1-NSP2, with a 99.6% sequence identity shared with CZE/H186/2018 and CZE/H187/2018, and were therefore clustered into lineage-III(Suppl. Fig. S1) . This result indicates that Czech G9P[4] strains contain Wa-like segment in the N1-NSP2 region, a trait reminiscent of G9P[4] strains found in Indonesia. To elucidate the origin of N1-NSP2, we collected the N1 sequences available from the GenBank database (Suppl. Fig. S1) .
The remaining 8 phylogenetic trees are shown in supplementary Fig. S3-S10 associated with the VP1-VP3, VP6, and NSP1-NSP5 genes. Indonesian reassortant G9P[4]/P[6]-DS-1-like strains exhibited the most significant genetic closeness to the G9P[4] DS-1-like strains detected in the Czech Republic in 2018, along with other contemporary DS-1-like RVA strains. Notably, these Indonesian strains differed slightly from G9P[4]-DS-1-like strains identified during the period from 2011 to 2015. Furthermore, the VP1-VP3, VP6, and NSP1-NSP5 genes in most of the Indonesian reassortant G9 DS-1-like strains showed genetic variations different from those observed in the equine-like G3P[8] strains found in Indonesia during 2015-2016. Notably, three strains detected in 2021-2022 (SOEP-733, SOEP-735, and SOEP-740) in this study showed similar genetic distances in the NSP3 and NSP4 genes to the equine-like G3P[8] strains.
Bayesian evolutionary analysis
To understand the origin and emergence timeline for unusual G9P[4] strains in Indonesia, we performed Bayesian evolutionary analysis on the G9-VP7, P[4]-VP4 and N1-NSP2 genes (Fig.3, Fig. 4, Suppl. Fig. S2). This analysis of the sequence data over time from the GenBank database enabled the estimation of substitution rates and the time of tMRCA for the G9P[4] genotype. Maximum clade credibility (MCC) trees were constructed using the Bayesian MCMC framework. The estimated time of tMRCA for the VP7 MCC tree including G9 was in 2004 with a 95% highest probability density (HPD) interval of 2001-2021.The estimated evolution rate was 6.43 × 10-3(4.12-8.74 × 10-3) (nt substitutions/site/year)(Fig. 3). The tMRCA for the P[4]-VP4 strains was 1994 (95% HPD interval 1990-2021) and the estimated evolution rate was 2.53 × 10-3 (1.54-3.52 × 10-3) (Fig. 4). The tMRCA of the P[6]-VP4 strains was 2008 (95% HPD interval 2000-2021), and the estimated evolution rate was 2.11 × 10-3 (1.35-2.86 × 10-3) (Fig. 4). Within the G9P[4] strains, the NSP2 gene of the reassortant G9-N1 strains diverged from the typical DS-1-like G9-N2 strains in 2003 (95% HPD interval 2000-2021) with the estimated evolution rate of 7.05 × 10-3 (4.22-9.88 × 10-3) and in 1998 (95% HPD interval 1996-2021) with the estimated evolution rate of 1.85 × 10-3 (1.48-2.22 × 10-3) (Suppl. Fig. S2). Based on the NSP2 MCC tree, the G9-N1 reassortment started in 2011 and was found in several countries such as the Czech Republic in 2018 and in India in 2013. The 11 segments of the unusual G9P[4] Indonesian strain had high nucleotide identities with all the segments from the Czech Republic (VP1-4, VP6, VP7, NSP1-5) ranging from 99.6-99.8% and were clustered into the same lineage.
Amino acid (aa) substitution of VP7 and VP6 Genes