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