Whole genomic analysis of G2P[4] human Rotaviruses in Mymensingh, north-central Bangladesh.

Rotavirus A (RVA) is a dominant causative agent of acute gastroenteritis in children worldwide. G2P[4] is one of the most common genotypes among human rotavirus (HRV) strains, and has been persistently prevalent in South Asia including Bangladesh. In the present study, whole genome sequences of a total of 16 G2P[4] HRV strains (8 strains each in 2010 and 2013) detected in Mymensingh, north-central Bangladesh were determined. These strains had typical DS-1-like genotype constellation. Most of gene segments from DS-1 genogroup exhibited high level sequence identities to each other (>98%), while slight diversity was observed for VP1, VP3, and NSP4 genes. By phylogenetic analysis, individual RNA segments were classified into one (V) or two-three lineages (V-VI or V-VII). In terms of lineages (sublineages) of 11 gene segments, the 16 Bangladeshi strains could be further classified into four clades (A-D) containing 8 lineage constellations, revealing the presence of three clades (A-C) with three lineage constellations in 2010, and a single clade (D) with four constellations in 2013. Therefore, co-existence of multiple G2P[4] HRV strains with different lineage constellations, and change in clades for the study period were demonstrated. Although amino acids in the antigenic regions on VP7 and VP4 were mostly identical to those of global G2P[4] strains after 2000, VP4 of clade D RVAs in 2013 had alanine and proline at positions 88 and 114, respectively, which are novel substitutions compared with recent global G2P[4] strains. Replacement of lineage constellations associated with unique amino acid changes in the antigenic region in VP4 suggested continuous genetic evolutionary state for emerging new G2P[4] rotavirus strains in Bangladesh.

Introduction

Rotavirus A (Group A rotavirus, RVA) is the leading etiological agent of severe gastroenteritis in infants and young children worldwide, and is estimated to cause 197,000 deaths in children <5 years of age (Lanata et al., 2013). As enteric pathogens, RVAs circulate in mammals and birds. Rotavirus is a genus of the family Reoviridae, and its genome is composed of 11 segments of double-stranded RNA enclosed in a triple-layered capsid. RNA segments of RVA encode six structural proteins (VP1-VP4, VP6 and VP7) and six nonstructural proteins (NSP1-NSP6) (Estes and Greenberg, 2013). Due to the segmented nature of the rotavirus genome, reassortment is considered to occur occasionally when co-infection with more than one rotavirus genotype occurs in hosts (Greenberg et al., 1981; Midthun et al., 1987; Urasawa et al., 1986). This reassortment has been revealed by whole genomic analysis (Ghosh and Kobayashi, 2011a).

The outermost layer of the rotavirus particle consists of two structural proteins VP7 and VP4, on which neutralization antigens are present and define independent serotypes (G and P serotypes, respectively). Based on the VP7 and VP4 genes, RVA has been genetically classified into G type and P type, respectively (Estes and Greenberg, 2013). A total of 27 G types and 37 P types have been described to date for human and animal rotaviruses (Matthijnssens et al., 2011; Trojnar et al., 2013). Genotypes G1P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8], G2P[4], G3P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8], G4P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8], G9P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8], G12P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8] and combinations thereof, are frequently detected in human RVAs throughout the world (Santos and Hoshino, 2005; Dóró et al., 2014). A comprehensive genetic analysis of all 11 segments revealed that there are two major genotype constellations in human rotaviruses (HRVs). The Wa genogroup (Wa-like genotype constellation) includes strains with G1P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8], G3P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8], or G4P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8] genotypes, and the DS-1 genogroup (DS-1-like genotype constellation) is usually associated with G2P[4] HRV strains (Matthijnssens et al., 2008).

G2P[4], one of the most common HRV genotypes worldwide, has been showing relatively high detection rates (16–36%) in South Asia recently (Dóró et al., 2014; Miles et al., 2012; Mullick et al., 2014). In Bangladesh, G2 has occasionally been the most common genotype among circulating HRVs for the past 20 years, while G9 and G12 emerged as predominant strains in the last decade (Afrad et al., 2013; Afrad et al., 2014; Ahmed et al., 2012; Dey et al., 2009; Miles et al., 2012; Paul et al., 2008). The reason for the occasional dominance of G2P[4] HRV in Bangladesh has not yet been fully understood, although the occurrence of multiple lineages of VP7 gene is suggested to be one of the possible reasons (Afrad et al., 2014).

G2P[4] genotype has been known to increase or persist in dominance after introduction of monovalent rotavirus vaccine in Brazil, Argentina, Australia, Belgium, and Korea (Gurgel et al., 2014; Kim et al., 2014; Kirkwood et al., 2009; Kirkwood et al., 2011; Mandile et al., 2014; Matthijnssens et al., 2014). Therefore, it is important to understand diversity and genetic evolution of G2P[4] HRVs in nature and their association with efficacy of vaccines. Whole genomic analysis of HRVs revealed presence of three distinct clades among G2P[4] in 2010–2011 winter season in the USA (Dennis et al., 2014). Furthermore, long-term investigations of whole genome of G2P[4] HRVs in Italy (Giammanco et al., 2014) and Japan (Doan et al., 2015) indicated occurrence of major genomic change in the global G2P[4] RVAs in the early 2000s.

In Bangladesh, genetic diversity and evolution of G2P[4] RVAs was analyzed in terms of VP7 gene from a large number of G2P[4] strains (Afrad et al., 2014). In this study, Afrad and coworkers revealed that multiple lineages of G1 and G2 HRVs had been co-circulating over the years. Whole genomes of HRV G1, G2 and G12 strains in Bangladesh were previously analyzed (Ghosh et al., 2011b; Rahman et al., 2007; Rahman et al., 2010), with only two G2P[4] HRV strains available for study (Ghosh et al., 2011b). Thus, the genomic diversity of G2P[4] HRVs in Bangladesh have not yet been well characterized using whole genome sequencing. The purpose of the present study was to elucidate evolutionary state of all the gene segments of G2P[4] HRVs, i.e., to understand genetic diversity of individual genome segments and correlation of evolution among gene segments. In the present study, we analyzed whole genome of 16 G2P[4] HRV strains in 2010 and 2013 to obtain clues to understand their persistence in Bangladesh.

Material and methods

Virus strains

A total of 17 G2P[4] HRV strains isolated from diarrheal stool samples collected from patients aged 3 months to 24 years in Mymensingh, Bangladesh were studied. Nine and eight HRV strains were obtained in Jan.–Feb. 2010 and Aug.–Dec. 2013, respectively (Table 1). In 2010, G2 was the most prevalent (41%), followed by G1 (25%), and G9 (8%), in Mymensingh (unpublished data). Genotyping results of HRV in 2013 have not yet been available. All the strains were confirmed to have G2P[4] genotypes by semi-nested RT-PCR as described previously (Iturriza-Gómara et al., 2004; Nagashima S et al., 2010). For the experiments, utmost attention was paid to avoid contamination especially in handling stool specimens, for example, using only one sample per day for RNA extraction and RT-PCR.

Table 1

Date, age, and sex of the patients infected with RVA strains.

RVA strain	Date of collection (year/month)	Age	Sex
RVA/Human-wt/BGN/J331/2010/G2P[4]*	2010.2	1Y	M
RVA/Human-wt/BGN/J306/2010/G2P[4]	2010.2	24Y	M
RVA/Human-wt/BGN/J303/2010/G2P[4]	2010.2	12Y	M
RVA/Human-wt/BGN/J300/2010/G2P[4]	2010.2	6Y	M
RVA/Human-wt/BGN/J266/2010/G2P[4]	2010.1	20Y	M
RVA/Human-wt/BGN/J265/2010/G2P[4]	2010.1	6M	M
RVA/Human-wt/BGN/J263/2010/G2P[4]	2010.1	20Y	M
RVA/Human-wt/BGN/J253/2010/G2P[4]	2010.1	5Y	M
RVA/Human-wt/BGN/J251/2010/G2P[4]	2010.1	22Y	M
RVA/Human-wt/BGN/M334/2013/G2P[4]	2013.9	8M	F
RVA/Human-wt/BGN/M315/2013/G2P[4]	2013.8	4M	M
RVA/Human-wt/BGN/M313/2013/G2P[4]	2013.8	10M	M
RVA/Human-wt/BGN/M312/2013/G2P[4]	2013.8	2Y	M
RVA/Human-wt/BGN/M310/2013/G2P[4]	2013.8	3Y	M
RVA/Human-wt/BGN/M292/2013/G2P[4]	2013.12	3M	M
RVA/Human-wt/BGN/M289/2013/G2P[4]	2013.11	5M	M
RVA/Human-wt/BGN/M282/2013/G2P[4]	2013.11	6M	F

*This strain had Wa-like genotypes in VP3 gene (M1) and NSP5/6 gene (H1) in the DS-1-like genotype constellation. Because the possibility of mixed infection could not be excluded, this strain was not included in the phylogenetic analysis.

Nucleotide sequencing, genotyping and sequence analyses

Viral RNA was extracted from stool samples using the QIAamp Viral RNA Mini Kit (Qiagen Science, MD, USA). RT-PCRs were performed using Superscript II reverse transcriptase (Invitrogen, Carlsbad, CA) and Prime Star GXL DNA polymerase (TaKaRa, Japan). Primers used for amplification of viral genome segments have been described previously (Ghosh et al., 2010a; Ghosh et al., 2010b; Ghosh et al., 2011b). For the RT-PCR of VP1-4, NSP1, and NSP4 genes, additional primers were designed in the present study as shown in Table 2. Nucleotide sequences were determined by the Sanger method using the BigDye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems, Foster City, CA) on an automated DNA sequencer (ABI PRISM 3100). The Basic Local Alignment Search Tool (BLAST) (http://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to search for the most similar rotavirus gene sequence and assign its genotype based on cut-off values indicated by Rotavirus Classification Working Group (Matthijnssens et al., 2011). Phylogenetic trees of individual gene segments were constructed by Maximum Likelihood method using MEGA6 software (Tamura et al., 2013), with sequences selected from GenBank database. Phylogenetic trees were statistically supported by bootstrapping with 1000 replicates, and phylogenetic distances were measured by the Kimura two-parameter model (Kimura, 1980). Multiple alignments of sequences were performed using CLUSTAL W ver. 2.1 program available on website of DDBJ (http://clustalw.ddbj.nig.ac.jp/). Sequence identity of a pair of gene sequences was determined by using LALIGN program on web server (http://www.ch.embnet.org/software/LALIGN_form.html).

Table 2

Primers designed in the present study.

Gene segment	Primer	Sequence(5'–3')	NucleotidePosition*	Size
VP1	390R	CAT CAA TGA GTC AGT GTA TTC	408–388	21 mer
VP1	2640R	GGT TTT ATG TCT TTA AGT ATG TCG	2666–2643	24 mer
VP2	480F	GGG GAC TAT GAT GTG AGA GAG	485–505	21 mer
VP2	306F	CCA ACA TTC GAA CCT AAA GAG ACG	299–322	24 mer
VP2	1014F	GGC GAG ATC GGT AGT ACC AG	997–1016	20 mer
VP3	2195R	GTA CCA CAT CTC ACA TTT GGC G	2195–2216	22 mer
VP3	1860R	CAC ATG TCC AGA CAC TGA ATT CTC	1888–1865	24 mer
VP3	2427R	TCG TGA TTG TCC AAA CGT GAT G	2425–2404	22 mer
VP3	1725R	CCC ATA TGA TTT GCA TAT TGA TC	1752–1730	23 mer
VP3	2124F	ATA TAG TAT AAC TTA TGC TGA CG	2113–2135	23 mer
VP4	2091F	GGA TAC ACT TAA TGA GAT CCC	2091–2111	21 mer
VP4	1215F	CTA TTA TGA ATG GCG GTGCTG	1190–1210	21 mer
NSP1	200R	ATG TTG ACA ACA ATC TAA GC	175–156	20 mer
NSP1	600F	ATG TAT TAC TGC TAG ATA GC	576–595	20 mer
NSP4	550F	AGA GGT TGA GCT GCC GTC GTC	569–589	21 mer
NSP4	500R	TAG CGT TTT CAC GTT CTT TTG	509–489	21 mer

*These positions correspond to individual gene segments of G2P[4] HRV.

Nucleotide sequence accession numbers

The GenBank accession numbers for the nucleotide sequences determined in the present study were listed in Table 3.

Table 3

GeneBank accesseion numbers for the 11gene segments of 16 G2P[4] RVA strains in Mymensingh, Bangladesh.

Strain	Viral gene segment
	VP1	VP2	VP3	VP4	VP6	VP7	NSP1	NSP2	NSP3	NSP4	NSP5
RVA/Human-wt/BGN/J306/2010/G2P[4]	KU199270	KU199271	KU199272	KU199273	KU199274	KU199275	KU199276	KU199277	KU199278	KU199279	KU199280
RVA/Human-wt/BGN/J303/2010/G2P[4]	KU248372	KU248373	KU248374	KU248375	KU248376	KU248377	KU248378	KU248379	KU248380	KU248381	KU248382
RVA/Human-wt/BGN/J300/2010/G2P[4]	KU248383	KU248384	KU248385	KU248386	KU248387	KU248388	KU248389	KU248390	KU248391	KU248392	KU248393
RVA/Human-wt/BGN/J266/2010/G2P[4]	KU248394	KU248395	KU248396	KU248397	KU248398	KU248399	KU248400	KU248401	KU248402	KU248403	KU248404
RVA/Human-wt/BGN/J265/2010/G2P[4]	KU248405	KU248406	KU248407	KU248408	KU248409	KU248410	KU248411	KU248412	KU248413	KU248414	KU248415
RVA/Human-wt/BGN/J263/2010/G2P[4]	KU248416	KU248417	KU248418	KU248419	KU248420	KU248421	KU248422	KU248423	KU248424	KU248425	KU248426
RVA/Human-wt/BGN/J253/2010/G2P[4]	KU356574	KU356575	KU356576	KU356577	KU356578	KU356579	KU356580	KU356581	KU356582	KU356583	KU356584
RVA/Human-wt/BGN/J251/2010/G2P[4]	KU356585	KU356586	KU356587	KU356588	KU356589	KU356590	KU356591	KU356592	KU356593	KU356594	KU356595
RVA/Human-wt/BGN/M334/2013/G2P[4]	KU199281	KU199282	KU199283	KU199284	KU199285	KU199286	KU199287	KU199288	KU199289	KU199290	KU199291
RVA/Human-wt/BGN/M315/2013/G2P[4]	KU356596	KU356597	KU356598	KU356599	KU356600	KU356601	KU356602	KU356603	KU356604	KU356605	KU356606
RVA/Human-wt/BGN/M313/2013/G2P[4]	KU356607	KU356608	KU356609	KU356610	KU356611	KU356612	KU356613	KU356614	KU356615	KU356616	KU356617
RVA/Human-wt/BGN/M312/2013/G2P[4]	KU356618	KU356619	KU356620	KU356621	KU356622	KU356623	KU356624	KU356625	KU356626	KU356627	KU356628
RVA/Human-wt/BGN/M310/2013/G2P[4]	KU356629	KU356630	KU356631	KU356632	KU356633	KU356634	KU356635	KU356636	KU356637	KU356638	KU356639
RVA/Human-wt/BGN/M292/2013/G2P[4]	KU356640	KU356641	KU356642	KU356643	KU356644	KU356645	KU356646	KU356647	KU356648	KU356649	KU356650
RVA/Human-wt/BGN/M289/2013/G2P[4]	KU356651	KU356652	KU356653	KU356654	KU356655	KU356656	KU356657	KU356658	KU356659	KU356660	KU356661
RVA/Human-wt/BGN/M282/2013/G2P[4]	KU356662	KU356663	KU356664	KU356665	KU356666	KU356667	KU356668	KU356669	KU356670	KU356671	KU356672

Results

The genotype constellation of 16 strains in this study was G2-P[4]-I2-R2-C2-M2-A2-N2-T2-E2-H2, showing typical DS-1-like genotype constellation. One strain J331 detected in 2010 had the genotype constellation G2-P[4]-I2-R2-C2-M1-A2-N2-T2-E2-H1, which contained genotypes of Wa genogroup (M1 and H1). Because the possibility of mixed infection with Wa genogroup HRV could not be excluded, strain J331 was not included in the phylogenetic analysis with other G2P[4] strains.

VP2, VP4, VP6, VP7, NSP1-3 and NSP5 genes from the 16 G2P[4] strains exhibited high level sequence conservation with >98% sequence identity to each other. In contrast, slight sequence diversity was observed for VP1, VP3, and NSP4 genes among the 16 strains (94.4%, 90.6%, and 94.3% sequence identity, respectively). Phylogenetic analysis together with representative G2P[4] RVAs from the world (Bányai et al., 2011; Chaimongkol et al., 2012; Doan et al., 2015; Ghosh et al., 2011b; Giammanco et al., 2014; Page and Steele, 2004) indicated that most of the Bangladeshi G2 HRV genome segments were classified into a single lineage (lineage V), while VP3 and NSP4 genes were assigned to two (V, VI) and three lineages (V–VII), respectively (Fig. 1, Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8, Fig. 9, Fig. 10, Fig. 11). We followed designation of lineages in each RNA segment as described by Doan et al. (2015).

Fig. 1

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding VP1. Bangladeshi RVA strains detected in 2010 and 2013 analyzed in the present study are marked with closed triangles and circles, respectively. Lineages and sublineages within a lineage are shown with vertical lines on the right. Lineages I–IV of individual genes were assigned by the scheme described by Doan et al. (2015), while other lineages and sublineages were designated in the present study. Scale bars are shown below. Bootstrap values are indicated at nodes of branches. Bootstrap values less than 70% are not shown.

Fig. 2

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding VP2. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 3

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding VP3. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 4

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding VP4. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 5

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding VP6. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 6

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding VP7. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 7

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding NSP1. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 8

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding NSP2. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 9

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding NSP3. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 10

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding NSP4. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Fig. 11

Phylogenetic dendrograms based on full-length nucleotide sequences of genes encoding NSP5. See legends of Fig. 1 for marks, lineage assignment, scale bars and bootstrap values.

Within the lineage V, if the sequences of present Bangladeshi HRVs clustered in a branch supported with high bootstrap value, sublineage was designated with a subscript attached with lineage V (Fig. 1, Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8, Fig. 9, Fig. 10, Fig. 11). A subscript “a” was assigned to five strains in 2013 because they clustered together in all the gene segments, and subscript b and c were assigned arbitrarily. G2P[4] strains in 2013 had mostly gene segments belonging to sublineage Va, while several genes of three strains were not assigned to lineage Va. In contrast, G2P[4] strains in 2010 were more divergent than 2013 strains, because individual gene segments belonged to various sublineages or were not assigned to the sublineages. Based on combination of sublineages (lineages) in 11 gene segments, G2P[4] HRVs in 2010 and 2013 were classified into four lineage constellations each (Table 4). In terms of the similarity of lineage constellations (five or more identical sublineages of lineage V, and lineages VI and VII), the eight strains each in 2010 and 2013 were classified into three clades (A-C) and one clade (D), respectively.

Table 4

Lineages (sublineages) of 11 genome segments detected in Mymensingh, Bangladesh.

RVA strain	Lineage (sublineages)* of viral protein genes (genotype : G2-P[4]-I2-R2-C2-M2-A2-N2-T2-E2-H2)											Lineage constellation#	Clade#
	VP7	VP4	VP6	VP1	VP2	VP3	NSP1	NSP2	NSP3	NSP4	NSP5/6
RVA/Human-wt/BGN/MMC6/2005/G2P[4]	V	Vc	Vb	Vb	Vc	VII	Vb	Vb	Vc	VI	Vb
RVA/Human-wt/BGN/MMC88/2005/G2P[4]	V	Vb	V	Vb	Vb	Va	V	V	V	VII	V
RVA/Human-wt/BGN/J306/2010/G2P[4]	V	Vb	V	Va	Vb	Va	V	V	Vb	V	V	2010–4	A
RVA/Human-wt/BGN/J303/2010/G2P[4]	V	Vb	V	Vb	Vb	Va	V	V	V	VII	V	2010–2	B
RVA/Human-wt/BGN/J300/2010/G2P[4]	V	Vb	V	Vb	Vb	Va	V	V	V	VII	V	2010–2	B
RVA/Human-wt/BGN/J266/2010/G2P[4]	V	Vb	V	Va	Vb	Va	V	V	Vb	Va	V	2010–1	A
RVA/Human-wt/BGN/J265/2010/G2P[4]	V	Vb	V	Va	Vb	Va	V	V	Vb	Va	V	2010–1	A
RVA/Human-wt/BGN/J263/2010/G2P[4]	V	Vb	V	Va	Vb	Va	V	V	Vb	Va	V	2010–1	A
RVA/Human-wt/BGN/J253/2010/G2P[4]	V	Vc	Vb	Vb	Vc	VI	Vb	Vb	Vc	VI	Vb	2010–3	C
RVA/Human-wt/BGN/J251/2010/G2P[4]	V	Vc	Vb	Vb	Vc	VI	Vb	Vb	Vc	VI	Vb	2010–3	C
RVA/Human-wt/BGN/M334/2013/G2P[4]	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	2013–1	D
RVA/Human-wt/BGN/M315/2013/G2P[4]	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	2013–1	D
RVA/Human-wt/BGN/M313/2013/G2P[4]	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	2013–1	D
RVA/Human-wt/BGN/M312/2013/G2P[4]	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	2013–1	D
RVA/Human-wt/BGN/M310/2013/G2P[4]	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	Va	2013–1	D
RVA/Human-wt/BGN/M292/2013/G2P[4]	Va	V	Va	V	Vb	Va	Va	Va	V	Va	Va	2013–2	D
RVA/Human-wt/BGN/M289/2013/G2P[4]	Va	Va	Va	V	Va	Va	Va	Va	Va	Vb	Va	2013–3	D
RVA/Human-wt/BGN/M282/2013/G2P[4]	V	Va	Va	V	Va	Va	Va	Va	V	Vb	Va	2013–4	D

Lineage (sublineage) is the same as that described in Fig. 1.

Lineage constellation was designated based on combination of lineages/sublineages of 11 genome segments, and clades A-D were defined as group of lineage constellations having five or more identical sublineages of lineage V (and lineage VI and VII in some strains).

All the genome segments of HRV strains analyzed in the present study clustered with those of G2P[4] HRV strains MMC6 and/or MMC88 detected in Bangladesh in 2005 (Ghosh et al., 2011b) and other contemporary G2P[4] strains from Americas, Asia, etc., within the lineage V, in the phylogenetic trees. However, clustering patterns with strains MMC6 and/or MMC88 were different depending on gene segments. Strains MMC6 and MMC88 had VP1, VP2, VP4, VP6, NSP1, NSP2, NSP5/6 genes classified into sublineages Vb or Vc. However, only VP3 gene of MMC88 was assigned to sublineage Va. VP7 gene of MMC6 and MMC88 did not cluster with the lineage Va. NSP4 genes of lineage VI (MMC6) and VII (MMC88) clustered with clade C and B strains, respectively.

Between sublineages Va and Vb, sequence identity of VP1 genes was approximately 94%, while >97.9 identity was found within the same sublineages. Lineage VI and VII NSP4 genes exhibited 94% identity to sublineage Va NSP4 genes. VP3 genes of sublineage Va showed 90.6% identity to those of lineage VI strains, in contrast to >99% identity among Va lineage.

Amino acid residues in VP7 defining neutralization domain of the G2P[4] strains in 2010 and 2013 were mostly identical to those of MMC88 strain in 2005, although one or two amino acid difference was found with strains J251 and J253 in 2010, and M313 and M289 in 2013 (Table 5). All the G2P[4] strains analyzed had 4–7 amino acids in the neutralization domain which are different from those of G2 component of pentavalent vaccine (RotaTeq) and G2 prototype strain DS-1, as described previously for recent G2P[4] rotaviruses (Afrad et al., 2014; Giammanco et al., 2014; Zeller et al., 2011). Although VP4 neutralization domains of G2P[4] HRVs in 2010 and 2013 were also similar to those of MMC88 strain, all the eight strains in 2013 had different amino acids at position 88 (A) and 114 (P) from those in the G2P[4] strains in 2005 and 2010 as well as DS-1 strain (Q and T, respectively) (Table 6a, Table 6b).

Table 5

Alignment of the amino acid residues defining the neutralization domains of VP7 (7-1a, 7-1b, and 7–2) between RotaTeq™ G2 component and G2 strains DS-1, TB-Chen, MMC88, and 16 G2P[4] RVA strains in Bangladesh.

RVA strain	7-1a															7-1b							7-2
	87	91	94	96	97	98	99	100	104	123	125	129	130	291		201	211	212	213	238	242		143	145	146	147	148	190	217	221	264
RVA/Vaccine/USA/RotaTeq SC2-9/1992/G2P[5]	A	N	S	D	E	W	E	N	Q	D	T	M	N	K		Q	D	V	S	N	S		R	D	N	T	S	D	I	S	G
RVA/Human-tc/USA/DS-1/1976/G2P[4]	A	N	S	D	E	W	E	N	Q	D	T	M	N	K		Q	D	V	D*	N	S		R	D	N	T	S	D	I	S	G
RVA/Human-tc/CHN/TB-Chen/1996/G2P[4]	A	N	S	D	E	W	E	N	Q	D	N	V	N	K		Q	D	V	N*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/MMC88/2005/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J306/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J303/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J300/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J266/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J265/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J263/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/J253/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	N	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	T*	G
RVA/Human-wt/BGN/J251/2010/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	N	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M334/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M315/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M313/2013/G2P[4]	T*	I*	S	N*	V*	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M312/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M310/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M292/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G
RVA/Human-wt/BGN/M289/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	V*	A*	G
RVA/Human-wt/BGN/M282/2013/G2P[4]	T*	N	S	N*	E	W	E	N	Q	D	T	M	D*	K		Q	D	V	D*	N	N*		R	D	N	T	S	D	I	S	G

Residues that differ from those of RotaTeq™ (G2 component).

Table 6a

Alignment of the amino acid residues defining the neutralization domains of VP4, VP8 subunit (8–1, 8–2, 8–3 and 8–4) between the P[4] strains (DS-1, TB-Chen, MMC88) and 16 G2P[4] RVA strains in Bangladesh.

RVA strain	8-1												8-2			8-3										8-4
	100	146	148	150	188	190	192	193	194	195	196		180	183		113	114	115	116	125	131	132	133	135		87	88	89
RVA/Human-tc/USA/DS-1/1976/G2P[4]	D	S	H	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	N	D		N	T	D
RVA/Human-tc/CHN/TB-Chen/1996/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	N	D		N	T	N
RVA/Human-wt/BGN/MMC88/2005/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J306/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J303/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J300/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J266/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J265/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J263/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J253/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/J251/2010/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	Q	T	N	N	E	N	S*	D		N	T	D
RVA/Human-wt/BGN/M334/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M315/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M313/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M312/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M310/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M292/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M289/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D
RVA/Human-wt/BGN/M282/2013/G2P[4]	D	S	Q	D	S	T	D	L	N	N	I		T	A		S	P*	T	N	N	E	N	S*	D		N	A*	D

Residues that differ from those of strains DS-1 and TB-Chen.

Table 6b

Alignment of the amino acid residues defining the neutralization domains of VP4, VP5 subunit (5–1, 5–2, 5–3, 5–4 and 5–5) between the P[4] strains (DS-1, TB-Chen, MMC88) and 16 G2P[4] RVA strains in Bangladesh.

RVA strain	5-1									5-2		5-3		5-4		5-5
	384	386	388	393	394	398	440	441		434		459		429		306
RVA/Human-tc/USA/DS-1/1976/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-tc/CHN/TB-Chen/1996/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/MMC88/2005/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J306/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J303/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J300/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J266/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J265/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J263/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J253/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/J251/2010/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M334/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M315/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M313/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M312/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M310/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M292/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M289/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R
RVA/Human-wt/BGN/M282/2013/G2P[4]	Y	F	L	W	P	G	R	T		P		T		L		R

Discussion

Recent molecular epidemiological studies of whole genome of G2P[4] HRVs for a long period conducted in Italy (Giammanco et al., 2014), and Japan (Doan et al., 2015) indicated that G2P[4] strains distributed globally appeared to have undergone intragenotype reassortment, observed by change of lineages of all the viral protein genes from the 1970s until 2011. As a result, it was suggested that major changes of genomic composition of the G2P[4] HRVs might occur in the early 2000s, when the “new” global G2P[4] HRV strains have been widespread replacing the “old” G2P[4] strains. In the present study, all the 16 RVA strains in 2010 and 2013 in Bangladesh belonged to the group of “new” G2P[4] HRVs clustering with recent global G2P[4] strains. This view was supported by the fact that amino acid residues in the VP7 antigenic regions were the same as those found in recent global G2P[4] HRVs, distinguishing them from old strains as well as the G2 component of pentavalent vaccine (Afrad et al., 2014; Dennis et al., 2014; Doan et al., 2011; Donato et al., 2014; Gómez et al., 2014; Zeller et al., 2011).

It was noted in the present study that change and replacement of sublineages and lineage constellations were observed for a short period of 3 years, although all the gene segments exhibited high sequence identities to each other. Individual viral genome segments were differentiated into 1–3 sublineages in a single lineage of individual genotypes (2 lineages in M2-VP3 and 3 lineages in E2-NSP4 genes). Thereby the 16 Bangladeshi strains were classified into four clades containing 8 lineage constellations, revealing the presence of three clades with four lineage constellations in 2010, and one clade with four constellations in 2013. Therefore, co-existence of multiple G2P[4] HRV strains with different lineage constellations was demonstrated in Bangladesh in the present study. Using classification of lineage (not sublineage as in the present study) of rotavirus genes, different allele constellations have been identified for G2P[4] HRV in a single winter season in a US community (Dennis et al., 2014), and also for G1, G3, and G4 RVAs in different study settings worldwide (McDonald et al., 2009; McDonald et al., 2011; McDonald et al., 2012; Wang et al., 2014).

Through analysis of VP7 from a number of G2 strains in the last 39 years in the US, lineage turnover was presumed to occur every 7 years on an average, generating new dominant strains (Dennis et al., 2014). Afrad and coworkers revealed that multiple lineages of G1 and G2HRVs co-circulated for one or a few seasons, followed by frequent replacement with different lineages, by the analysis of the VP7 gene of HRVs detected for more than 20 years in Bangladesh (Afrad et al., 2014). Therefore, lineage of the G2-VP7 gene is considered to change frequently in global level, as suggested by another study in Australia (Donato et al., 2014). Our present study, despite a short period, change of sublineages was documented for a whole genome, providing more detailed evidences than lineage.

In our study in Bangladesh, lineage groups (clade A-C) in 2010 were considered to be replaced by clade D in 2013. However, clade A had VP1, VP3, and NSP4 genes belonging to sublineage Va which was commonly found in clade D. Therefore, it is suggested that clade A HRVs in 2010 were one of the ancestral viruses to generate clade D HRVs in 2013. Viral genome segments of clade D HRV were mostly assigned to sublineage Va, which was unique to present Bangladeshi strains in 2013 without clustering with other global strains. This finding suggested that strains of clade D may be newly emerging G2P[4] HRVs in Bangladesh. Furthermore, VP4 of clade D HRVs characteristically possess alanine and proline at position 88 and 114 in the antigenic region, respectively. Both amino acids are not found in Bangladeshi strain MMC88 as well as global G2P[4] strains (Giammanco et al., 2014). Proline at position 114 is conserved in P[4] rotaviruses in Nepal in 2004-2005: a global and evolutionary perspective. Arch. Virol. 2011 ">8] VP4, and also identified in only a few G2P[4] Brazilian strains in 2010 (Gómez et al., 2014). It is therefore important to survey global spread and distribution of the G2P[4] RVA genes belonging to the clade D, and their antigenic characteristics of VP4. In contrast, clade C strains that were detected in only 2010 had asparagine at position 130 of VP7, like strains DS-1 and TB-Chen, suggesting that these strains have a trait of old G2HRVs.

VP3 gene of a Bangladeshi strain MMC88 in 2005 clustered with most G2P[4] HRV strains in the sublineage Va. The MMC88-VP3 gene is genetically related to caprine rotavirus strain GO34 from Bangladesh and suggests an animal origin (Ghosh et al., 2010a; Ghosh et al., 2011b), while another G2P[4] strain, MMC6, detected in 2005 has a VP3 gene frequently detected in HRVs. It was interesting in the present study that such caprine-like VP3 gene was identified in most of G2P[4] strains in 2010 and 2013 (lineage Va), suggesting successful adaptation of this gene to Bangladeshi HRVs of DS-1 genogroup. This may be evidenced by the fact that the VP3 genes clustering with the caprine strain and MMC88 have been detected in USA, Brazil, Thailand, Australia, and Italy (Dennis et al., 2014; Giammanco et al., 2014; Gómez et al., 2014).

In conclusion, our present study elucidated that multiple, genetically distinct G2P[4] HRVs are circulating in north-central Bangladesh, by whole genome-based phylogenetic analysis. Replacement of genomic constellations, and amino acid change in the antigenic region in VP4 were observed even for a short period from 2010 to 2013. Further continuous surveillance of G2P[4] HRVs is necessary at a global level to understand their evolutionary state.

Declarations

Author contribution statement

Satoru Aida, Nobumichi Kobayashi: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Samsoon Nahar, Shyamal K. Paul, Muhammad A. Hossain, Muhammad R. Kabir, Santana R. Sarkar, Salma Ahmed, Souvik Ghosh, Meiji S. Aung: Performed the experiments.

Noriko Urushibara, Mitsuyo Kawaguchiya, Ayako Sumi: Contributed reagents, materials, analysis tools or data.

Funding statement

Nobumichi Kobayashi was supported by the Grant-in-Aid for Scientific Research (Grant no. 25305022) from the Japan Society for the Promotion of Science.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.