The first RVA vaccine was approved in the United States in 2006. The World Health Organization (WHO) currently recommends the use of four oral RVA vaccines: Rotarix® (GlaxoSmithKline Biologicals, Rixenstart, Belgium), RotaTeq® (Merck Sharp & Dohme LLC, Rahway, NJ), RotaSiil® (Serum Institute of India, Pune, India), and Rotavac® (Bharat Biotech, Telangana, India)8. Clinical trials and real-world usage of these vaccines have validated the safety and efficacy of all RVA vaccines. However, their effectiveness is lower in developing countries than in developed countries. In low-income countries lacking government RVA vaccination programs, diarrhea caused by RV infection has serious consequences, contributing to the diversity of RVA genotypes9. In Indonesia, two RVA vaccines, Rotarix and RotaTeq, have been available commercially since 2013. However, neither of these vaccines has been included in the universal immunization program thus far.
RVA belongs to the family Reoviridae . RVA is adouble-stranded RNA virus consisting of 11 segments of the genome encoding six structural proteins (VP1, VP2, VP3, VP4, VP6, and VP7 genes) and six non-structural proteins (NSP1, NSP2, NSP3, NSP4, and NSP5/NSP6)2,3. The outer layer proteins consist of VP7 glycoprotein (the G genotype) and VP4 protein (the P genotype). Up to June 2024, Rotavirus Classification Working Group has proposed at least 42G, 58P, 32I, 28R, 24C, 23M, 39A, 28N, 28T, 32E, and 28H genotypes (https://rega.kuleuven.be/cev/viralmetagenomics/virus-classification/rcwg). The common recombinant genotypes were found in seven G genotypes (G1-G4, G8, G9 and G12) and three P genotypes (P[4], P[6], and P[8]), all of which are usually associated with human infections2,10,11. RVA genomic constellations are divided into three groups, Wa-like constellation (G1-P[8]-I1-R1-C1-M1-A1-N1-T1-E1-H1), DS-1-like (G2-P[4]-I2-R2-C2-M2-A2-N2-T2-E2-H2), AU-1-like (G3-P[9]-I3-R3-C3-M3-A3-N3-T3-E3-H3) a minor genomic constellation2,4,11.
The existence of the inter-genogroup reassortants highlights the wide diversity within RVA13. Moreover, recent findings on these inter-genogroup reassortant strains underscore the ongoing spread of rare RVA strains across Asia and other regions12. Although the reassortant RVA strains have been described in many studies, they appeared only sporadically until the emergence of reassortant G1P[8]-DS-1-like reported in Japan in 201413. We previously reported the equine-like G3 RVAs as predominant strains among children in Indonesia in 2015-2016. However, from July 2017 until 2018, the predominant RVA strains were completely replaced by typical human G1 and G3 genotypes, suggesting the dynamic changes in RVA genotypes from equine-like G3 to typical human G1 and G312.
The G9 genotype have emerged and spread globally since 1995. It is now acknowledged as the fifth major human RVA genotype14. The geographical distribution of G9 strains has contributed to their diversity, leading to the identification of various combinations involving the G9 genotype and the P genotype, including G9P[8], G9P[6], G9P[4], and G9P[11]13,15-17. It is noteworthy that the majority of G9 RVA reported cases were characterized as a combination of P[8] genotypes16. Reassortant G9P[4]-DS-1 like RVA strains exhibited notably high prevalence rates in several countries, including Mexico (80%), Guatemala (66%), Bangladesh (31.8%), India (35.3%), Iran (30.3%), and Ghana (16%)15,17-19. In Japan, G9P[4] prevalence has surged, reaching the levels between 11.5% and 21.2% during the years 2010-201213. Although unusual rotavirus strains have been described, they appeared only sporadically until the reassortment G1P[8] genotype on the DS-1 genetic backbone was discovered in several countries20. Notably, the emergence of G9 RVA with the P[8] genotype combination has been documented in many studies4,6,21.
During RVA surveillance conducted in Sidoarjo, Indonesia, we initially identified four cases of G9P[4] RVA infections in the years 2018 and 2020. Surprisingly, in 2021, the majority of RVA-positive specimens (19 out of 21) were identified as G9P[4], while the remaining three strains were G9P[6]. To gain insights into the genetic characteristics of these G9P[4]/P[6] strains, we conducted whole-genome sequencing for all 26 of the detected G9P[4]/P[6] strains using next-generation sequencing (NGS). In addition, we investigated the evolutionary patterns and genetic diversity among Indonesian G9P[4] strains, along with comparisons to RVA strains detected in other countries. We demonstrate that the emergence of G9P[4] was derived from the genetic diversity circulating in Indonesia, whereas the predominant strain survived for about a year and changed rapidly alongside or in connection with other genotypes.
MATERIAL AND METHODS
Specimens
We collected 356 stool samples from children under 5 years old who were hospitalized for acute gastroenteritis (AGE), in Sidoarjo, Indonesia between August 2018 and February 2022. AGE was defined as three times or more occurrence of looser than normal stools within 24 h. All stool samples were screened by an immunochromatography, Dipstick “Eiken” Rota kit (Eiken Chemical, Tokyo). The kit detected a total of 26 RVA-positive samples, which were analyzed by RT-PCR to determine their G and P genotypes. We primarily focused on detecting and characterizing G9P[4]/P[6] strains, which were subsequently subjected to whole-genome sequencing at the National Institute of Infectious Diseases in Tokyo.
This study was approved by the Research Ethics Board of both hospitals, Universitas Airlangga (ethics approval number 2054/ UN3.14/LT/2015) in Indonesia and Kobe University (ethics approval number 1857) in Japan. Written informed consent was obtained from the children’s parents or guardians.
RNA extraction and RT-PCR genotyping
To prepare a 10% (wt/vol) stool suspension, we utilized phosphate-buffered saline (PBS; pH 7.4) and clarified it by centrifugation at 21,130×g for 10 min. The supernatant was collected and stored at –80°C until further use. The viral RNA was extracted using the Direct-zol™-96 MagBead RNA (with DNase treatment) kit (Zymo Research, Irvine, CA) . In brief, 100 μl of stool suspension in PBS was mixed with 300 μl of Trizol LS reagent (Thermo Fisher Scientific, Waltham, MA). Further RNA extraction steps were then carried out using the automated KingFisher Flex System (Thermo Fisher Scientific). The extracted viral RNAs were eluted using DNase/RNase‐free water and were employed for RT‐PCR genotyping and whole-genome analysis by NGS.
All of the positive samples were subjected to genotyping in the VP7 (G typing) and VP4 (P typing) genes by multiplex RT-PCR. The primer sets of VP7 and VP4 were described previously 12, which enabled us to genotype the equine-like G3 strain along with other strains (G1, G2, typical human G3, G4, G8, G9, and G12)13. The RT‐PCR products were subjected to 2% agarose gel electrophoresis.
Complementary DNA (cDNA) library building and Illumina MiSeq sequencing
The cDNA library preparation and the Illumina Miseq sequencing were conducted as previously described15. Briefly, a 300‐bp fragment library ligated with bar‐coded adapters was constructed for individual strains using an NEBNext Ultra RNA library Prep Kit for Illumina ver. 1.2 (New England Biolabs, Ipswich, MA). Library purification was performed with Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA). The quality of the final DNA libraries was assessed using the High Sensitivity D1000 ScreenTape assay on 4150 TapeStation (Agilent Technologies, Santa Clara, CA) and quantified with Qubit™ 1X dsDNA High Sensitivity (HS) using the Qubit Flex Fluorometer (Thermo Fisher Scientific). Sequencing was performed on an Illumina MiSeq sequencer (Illumina, San Diego, CA) using the Miseq Reagent Kit ver. 2 (Illumina) to generate 151 paired end reads. Data analysis was performed using QIAGEN CLC Genomics Workbench v22 (CLC Bio, Tokyo, Japan). The nucleotide sequences of each gene segment of 26 Indonesian G9P[4]/P[6] strains was obtained by de novo assembly.
Phylogenetic analysis
Nucleotide sequence comparisons were carried out with the references retrieved from GenBank by using the BLAST program (http://blast.ncbi.nlm.nih.gov/) with the sequences of our G9P[4] and G9P[6] strains as the query sequences for the VP7 (G9 genotype), VP4 (P[4] and P[6] genotype) and NSP2 (N1 and N2 genotype) genome segments. Phylogenetic trees were constructed using Molecular Evolutionary Genetic Analysis (MEGA) ver. 11.0 with reference sequences retrieved from the DDBJ/EMBL/GenBank databases using the Maximum Likelihood method. Alignments were performed using the CLUSTAL W in MEGA 11 and phylogenetic trees were constructed by the neighbor-joining method. To confirm the reliability of phylogenetic tree analysis, bootstrap resampling and reconstruction were carried out 1,000 times22. The full-length genome sequences determined in this study were subjected to multiple alignment with the sequences obtained from the GenBank database using the MAFFT multiple sequence alignment, ver. 7.023.
Bayesian evolutionary analysis using BEAST
Estimation time of Most Recent Common Ancestor (tMRCA) was determined for the VP7, VP4 and NSP2 genes of the G9P[4] and G9P[6] strains by the Bayesian Markov chain Monte Carlo (MCMC) method in BEAST ver.1.8.124. The models used for BEAST analyses of G9P[4] and G9P[6] were GTR + G (VP7), T92 + 1 (VP4) and GTR + G (NSP2). Strict clock, relaxed clock and coalescent exponential growth models were used. MCMC runs were carried out for 100 million generations and evaluated using Tracer ver.1.624. Only parameters with an effective sample size of >200 were accepted. The maximum clade credibility (MCC) trees were annotated with the Tree annotator and viewed with FigTree ver.1.424.
Nucleotide sequence accession numbers
All the nucleotide sequences of the 11 gene segments of G9P[4] and G9P[6] have been deposited in the DDJB/GenBank/EMBL database under the accession numbers PP948908 - PP949193.
RESULTS
Nucleotide sequence and genotype constellation of Indonesian G9 strains