Molecular epidemiology and surveillance of circulating rotavirus among children with gastroenteritis in Bangladesh during 2014–2019

Acute gastroenteritis is one of the major health problems in children aged <5 years around the world. Rotavirus A (RVA) is an important pathogen of acute gastroenteritis. The burden of rotavirus disease in the pediatric population is still high in Bangladesh. This study investigated the prevalence of group A, B, and C rotavirus (RAV, RBV, RCV), norovirus, adenovirus (AdV) and human bocavirus (HBoV) infections in children with acute gastroenteritis in Bangladesh from February 2014 to January 2019. A total of 574 fecal specimens collected from children with diarrhea in Bangladesh during the period of February 2014-January 2019 were examined for RAV, RBV and RCV by reverse transcriptase- multiplex polymerase chain reaction (RT- multiplex PCR). RAV was further characterized to G-typing and P-typing by RT-multiplex PCR and sequencing method. It was found that 24.4% (140 of 574) fecal specimens were positive for RVA followed by AdV of 4.5%. RBV and RCV could not be detected in this study. Genotype G1P[8] was the most prevalent (43%), followed by G2P[4] (18%), and G9P[8] (3%). Among other genotypes, G9P[4] was most frequent (12%), followed by G1P[6] (11%), G9P[6] (3%), and G11P[25] (3%). We found that 7% RVA were nontypeable. Mutations at antigenic regions of the VP7 gene were detected in G1P[8] and G2P[4] strains. Incidence of rotavirus infection had the highest peak (58.6%) during November to February with diarrhea (90.7%) as the most common symptom. Children aged 4–11 months had the highest rotavirus infection percentage (37.9%). By providing baseline data, this study helps to assess efficacy of currently available RVA vaccine. This study revealed a high RVA detection rate, supporting health authorities in planning strategies such as introduction of RVA vaccine in national immunization program to reduce the disease burden.

Introduction All children aged <15 years with viral gastroenteritis (clinically suspected) were enrolled in the study. Samples were collected from pediatric patients with three or more watery nonbloody stools/24 hours. Non-infectious bloody diarrheal cases were excluded following WHO case standard guideline. The ages ranged from 55 days to 15 years with a mean age of 24 months. One fecal specimen was collected from each patient. Informed consent of the parents was taken before sample collection. Ten percent suspension of the fecal specimens were made and centrifuged at 10000x g for 10 min. All the specimens were stored at -20˚C until further analysis.

Ethical approval
This study was ethically approved by the Biosafety, Biosecurity & Ethical Committee (BBEC) of Jahangirnagar University. Ethics committee approval number for this study is BBEC, JU/M 2020/(10)3.

Viral RNA and DNA extraction
Fecal suspensions were thawed and centrifuged at 10000x g for 10 min. Viral RNA was extracted using 140 μL of the supernatant by SV Total RNA Isolation System (Promega kit) according to the manufacturer's protocols (Promega, Madison, USA). Viral DNA for adenoviruses and human bocaviruses was extracted following the established protocols in .

Complementary DNA (cDNA) synthesis by reverse transcription
For rotavirus and norovirus, cDNA was synthesized using Promega cDNA synthesis kit (Promega, Madison, USA). For reverse transcription (RT), 3 μL of extracted RNA (100 ng) was mixed with a reaction mixture, consisting of 1 μL of oligo dT primer (25 μg/ml) and 1 μL of nuclease free water in nanocentrifuge tube, kept at 70˚C for 5 min and chilled for 5 min in ice. After that 4 μL of 5X reaction buffer, 2 μL of 5 mM MgCl 2 , 1 μL of 10 mM PCR Nucleotide Mix, 0.5 μL of ribonuclease inhibitor (1 u/ μL), 1 μL of reverse transcriptase, and 6.5 μL of nuclease free water were mixed for each sample in the same nanocentrifuge tube, containing 5 μL of previously chilled genomic RNA. Total volume of the mixture was 20 μL. The mixture was heated at 25˚C for 5 min, 42˚C for 60 min and 70˚C for 15 min. The cDNA was prepared and stored at -20˚C.

Polymerase chain reaction (PCR)
For rotavirus A, two specific primers F-1 GGCTTTAAAAGAGAGAATTTC 28 and R-373 ACT GATCCTGTTGGCCATCCTTT 395 were used to amplify VP7 gene. For norovirus, two specific primers F-5342 CGCCGTGGCTCCTGCTCT 5361 and R-5671 GTTCGCCATCACAAAAGATGTG 5653 were used to amplify VP1 gene. For adenovirus, two specific primers AD1-F 1834 TTCCCCATGGCICAYAACAC 1853 and AD2-R 2315 CCCTGGTAKCCRATRTTGTA 2296 were used to amplify hexon gene. For human bocavirus, two specific primers AK-VP-F1 3274 CGCCGTGGCTCCTGCTCT 3291 and AK-VP-R1 3860 TGTTCGCCATCACAAAAGATGTG 3876 were used to amplify VP1 gene [25]. The PCR reaction was performed at 94˚C for 3 min, followed by 35 cycles of denaturation at 94˚C for 1 min, annealing at 55˚C for 1 min and extension at 72˚C for 1 min. Final extension was done at 72˚C for 7 min and then held at 4˚C. PCR reaction mixture contained 12.5 μL of the 2X master mix (GoTaq 1 Green Master Mix, Promega, USA), 1 μL of 1 μM forward (F) primer, 1 μL of 1 μM reverse (R) primer, 6.5 μL of nuclease free water, and 4 μL of template (200 ng). For rotavirus and norovirus, 4 μL cDNA (200 ng) was used as template. Total volume of the PCR reaction mixture was 25 μL. Previously sequenced samples were used as positive controls. The PCR was performed in 2720 Thermal Cycler (Applied Biosystems, USA) [25].

Group A rotavirus G typing
Group A rotavirus G typing was conducted using the protocol from the previously described method [23]. Reverse transcription (RT) was performed for the full-length VP7 gene, and the first amplification was completed by using Beg9 and End9 primers. The second amplification was performed using the first PCR product as the template with G type-specific mixed primers 9T1-1, 9T1-2, 9T-3P, 9T-4, and 9T-B and forward primer 9con1 (9con1 is specific primer that is used to characterize G genotype). Amplicons of 158 bp, 224 bp, 466 bp, 403 bp, and 110 bp were specifically generated for G1, G2, G3, G4, and G9, respectively.

Agarose gel electrophoresis
The amplicons were electrophoresed using 1.5% agarose gel. The gel run was continued for 25 min horizontally. Separated amplicons were visualized using UV spectrophotometer (SPE-CORD-205, Analytik-Jena, Germany) and photographed in the Gel doc system (BioRad, USA).

Nucleotide sequence analysis
The nucleotide sequences of PCR amplicons positive for rotavirus A, norovirus, human bocavirus and adenovirus were determined with the Big-Dye terminator cycle sequencing kit and an ABI Prism 310 Genetic Analyzer (Applied Biosystems Inc. Foster City, CA) [25]. The sequences were analyzed using Chromas 2.6.5 (Technelysium, Helensvale, Australia). The sequence data were converted into FASTA format. Sequence homology was confirmed using the BLASTn (https://blast.ncbi.nlm.nih.gov/Blast.cgi) program. Multiple sequence alignment was performed in the BioEdit 7.2.6 software using the ClustalW Multiple Alignment algorithm [26]. The similarity matrix was computed using the Maximum Composite Likelihood model. The nucleotide sequence data of all Bangladeshi RVA had been submitted to the GenBank database.

Phylogenetic analysis
Phylogenetic and evolutionary relationship analysis of RVA was conducted using VP7 gene and the reference sequences by the MEGA X software [27]. Phylogenetic trees with 1000 bootstrap replicates of the nucleotide alignment datasets were generated. For calculating the genetic distance, Kaimura-2 parameter model was used. Maximum Composite Likelihood (MCL) method was used for phylogenetic relationship analysis [28]. The reference strains used for RAV phylogenetic analysis are provided in the S1 Chart.

Nucleotide sequence submission
The nucleotide sequence data used in this study had been submitted to the GenBank database and had been assigned accession number MN202241 and MT310555-MT310574.

Statistical analysis
Statistical analysis was conducted using IBM 1 SPSS 1 Statistics for windows, version 23.0 (IBM Corp. Armonk, New York, USA) software package. Prevalence of diarrheal pathogens were detected by dividing number of cases by total number of study children during the study period. Statistical analysis for phylogenetic data was calculated by kimura-2-parameter model in the MEGA X software [27,28].

Coinfection of children with acute gastroenteritis
Frequencies of both monoinfection and coinfection of pediatric patients were determined. Infection of viruses was found in 75.4% (159 of 211) of the cases. Coinfection between viruses was found in 24.6% (52 of 211) of cases. Dual infection between RV and NoV accounted for most of the cases (10) of coinfection, followed by RV-AdV (6), RV-HBoV (4), NoV-AdV (3), NoV-HBoV (2) and AdV-HBoV (1).

Age and gender distribution of rotavirus among children with acute gastroenteritis
Viruses associated with acute gastroenteritis were detected in every group aged <60 months.

Clinical features of acute gastroenteritis associated with rotavirus infection
Viruses associated with acute gastroenteritis were found to be involved with various clinical complications of the pediatric population. Among rotavirus infected patients, diarrhea (90.7%) was the most common, followed by dehydration (88.7%), vomiting (84.3%), fever (77.8%) and abdominal pain (67.1%) (Fig 1). Whereas diarrhea (85.7%) was the most frequent in the norovirus infected children, vomiting (87.5%) was the most common in the adenovirus infected children, and dehydration (94.7%) was the most common in the human bocavirus infected children.

Seasonal pattern of rotavirus infection
Rotavirus infection was detected all year round. But majority of rotavirus infection cases (58.6%) were found during the winter season (Nov-Jan), followed by 30.7% in the rainy season (May-Jul), 7.1% in the spring season (Feb-Apr), and 3.6% in the autumn season (Aug-Oct), respectively. During the winter season the average temperature was around 17˚C. The incidence of rotavirus had a peak at a lower temperature in Bangladesh (Fig 2).

Distribution of G1 genotypes into unique sublineages
The sequences of VP7 genes of G1 isolates detected in this study and worldwide rotavirus reference G1 strains were used to investigate the heterogeneity of rotaviruses within G1 genotype. Twelve distinct lineages and 23 sublineages were identified (Fig 8).  (Fig 8). G1 isolates of lineage IV were closely related with G1 rotavirus of Finland and Russia. Within each sublineage, the nucleotide identities of rotavirus strains ranged from 97% to 100%, indicating less than 3% genetic difference among them. On the contrary, the nucleotide sequence variation of rotavirus G1 strains between lineages was higher, ranging from 6% to 12%.

Mutation analysis of VP7 gene of group A rotavirus
Rotavirus VP7 protein has three important antigenic regions (region A-amino acid position 87 to 100, region B-amino acid position 141-150, and region C-amino acid position 208-224), which can be variable within genotype. About 90%-99% sequence similarities were found in the genotype G1P [8] of this study with vaccine strain MF469224 and JN849114. In five G1P [8] strains, substitution point mutations were detected in the region A. However, all G2P [4] strains had the same mutation G216R at the region C comparing with vaccine strain GU565068. About 95%-99% sequence similarity was detected among G2P [4] strains. No significant mutation was found for G9P [4] strains in the regions A, B or C. Two G1P [6] strains had a deletion mutation at the position 220. G9P [8] strains had several mutations at regions A and C, sharing 95% sequence similarity. G9P [6] strains had substitution mutation at the regions B and C, sharing 96% sequence similarity. Isolate of G11P [25] RVA had only one mutation at the position T97I, sharing 100% sequence similarity with the reference (Table 3).

Still in this 21 st century rotavirus infection is a burden for children in developing countries like
Bangladesh. Virus infection in the gastrointestinal tract is one of the major causes of children morbidity and mortality in developing countries [1,2,4]. In agreement with previous studies, this study detected 36.8% (211 of 574) prevalence of virus infection in children with acute gastroenteritis in Bangladesh [19]. It was also found that rotavirus (RV) had the highest prevalence (24.4%), followed by norovirus (4.9%), adenovirus (4.2%) and human bocavirus (3.3%) [25]. Prevalence of viral infection in diarrheal children with AGE reported in this study were similar with the previous studies in developing countries like Nepal, India, Sri Lanka, Pakistan, South America, Africa and Europe [19][20][21][29][30][31][32][33]. In this study, rotavirus A was detected in 24.4% fecal specimens. No RVB or RVC was detected in this study from Bangladeshi children. This high prevalence of RVA in children with AGE in Bangladesh is an indication of lower vaccine efficacy or poor vaccine coverage. This finding of high RVA prevalence in children was consistent with the previous studies during pre-vaccinated era of rotavirus in Bangladesh and other developing countries [19][20][21][29][30][31][32][33]. Other possible causes of high prevalence of RVA might be poor hygienic conditions of most of the residents, presence of great genotypic diversity and not including vaccine in the national immunization program in Bangladesh [19]. In this study, the coinfection of other viruses with rotavirus was also determined. Coinfection between diarrheal viruses was detected in 24.6% of the positive cases. Coinfection between rotavirus and norovirus accounted for most of the cases (10 of 52) of coinfection. These findings are in good agreement with the previous studies of rotavirus coinfection in many developing countries worldwide [33][34][35]. Children aged <60 months were divided into six age groups according to previous study [25,33]. In this study, rotavirus was detected to be most prevalent (37.9%) among children   [20,[29][30][31][32][33][34][35]. Rotavirus infection has been reported to cause various clinical symptoms in diarrheal children worldwide. In this study, diarrhea (90.7%) dehydration (88.7%) and vomiting was found in higher percentage in rotavirus positive cases than negative cases. However, abdominal pain was reported in higher percentage from rotavirus negative cases. Clinical symptoms associated with rotavirus infection reported in this study also agreed with the previous studies in Bangladesh, Asia, Africa and Europe [20,21,[29][30][31][32][33].
While analyzing RVA infection seasonality, we found high incidence of rotavirus during the winter (Nov-Jan) season with low temperature in Bangladesh. About 58.6% of rotavirus infection in children were detected during the winter during 2014 to 2019. Another peak of rotavirus infection (30.7%) was found during the rainy season (May-Jul). During this study, the incidence of rotavirus increased with declining temperature and we found highest peak (21% infection) at 15˚C average temperature in Bangladesh. High prevalence and peak of RVA infection had been reported in Dey et al. (2009) and Phan et al. (2007). The seasonality of RVA in this study is with good agreement with studies in Bangladesh, South East Asia and Japan [20,36,37].
One of the major objectives of this study was to characterize the genotype of circulating rotavirus in Bangladesh. Changing pattern of both G and P genotypes along with their combinations were characterized. Worldwide epidemiological research on rotavirus specified that G1 is the predominant genotype. However, incidence of genotype G2-G4, G9 and G12 has also been detected over time [14,28]. During 1992 to 1997, rotavirus G4 genotype was the most common in Dhaka but became less common over time. During 2002-2005, genotype G1 was the most prevalent, while in 2005-2006, G2 was prevalent over G1 in Bangladesh [18][19][20]. Genotype G9 became predominant during 2006 to 2012 along with G1 in Bangladesh [19]. We found increased frequency of G1 genotype (57%) followed by G2 (20%), G9 (20%) and G11 (2.9%), respectively in children with acute gastroenteritis in Bangladesh. However, the frequency of G1 was found reducing in Giri et al. (2018) in developing countries of Asia. Interestingly, this study reported decreased frequency of G9 and increased frequency of G2 from previous studies, Mahmud-Al-Rafat et al. (2018) in Bangladesh [19]. Of note, G11 was also specified in this study. RVA studies, Giri et al. (2018) and Sadiq et al. (2019) have reported that prevalence of G3 has increased than before in developing countries in Asia and Europe in recent time [33]. However, genotype G3, G4 along with G12 was not detected in this study. Increased prevalence of G1 is in good agreement with previous reports but absence of G3, G4 and decreased prevalence of G9 along with increased incidence of G2 and G11 are new in Bangladesh [18][19][20][21]. This altered trends of RVA genotype distribution may also reflect regional variation in recent times. In this study, four P genotypes-P [8], P [4], P [6] and P [25] were detected. Genotype P [8] was the most frequent (45.7%) followed by P [4] (30%), P[6] (14.3%) and P [25] (2.9%) and 7.1% was nontypeable. Among P genotype, P [8] (80%) has remained the most prevalent genotype followed by P [4] and other non-P [8] genotype (20%) in previous studies in Bangladesh, India, Nepal and Pakistan [20,28]. Of note, a decrease of P [8] frequency was replaced by an increase frequency of non-P [8] genotypes in Bangladesh during 2014-2019. To the best of our knowledge, this will be the first report of genotype P [6] in Bangladesh.
There are no recent reports of P [6] from Bangladesh, but it has been reported from Nepal and Myanmar in Giri et al. (2018).
In the G and P genotype combination, G1P [8] [8] had been also reported in Bangladesh [19,35]. Another unusual genotype G9P [4] has emerged in Bangladesh in recent time and this genotype has also been reported in Asia, Europe, Africa and Americas [19,[29][30][31][32][33]. However, we also detected a great genotypic diversity of RVA during 2014-2019 in Bangladesh. Seven genotype combinations were revealed in this study. The most prevalent genotype was G1P [8] with the highest frequency of 43%, followed by G2P [4] (18%). This genotypic distribution of RVA in Bangladesh was in good agreement with the previous reports in India, Nepal, Pakistan, Japan, Brazil, Europe, Africa and Bangladesh [18-20, 23, 28-33]. Interestingly, we found that the frequency of previously predominant G9P [8] infection has reduced to 3%. We also detected high prevalence (12%) of genotype G9P [4], followed by genotype G1P [6] (11.4%), G9P [6] (3%), and G11P [25] (3%). To the best of our knowledge, it is the first report of G1P [6] and G9P [6] in Bangladesh. There are limited studies that had reported G1P [6] and G9P [6] in Nepal and Myanmar during 2009-2015 and Sadiq et al. (2019) has reported G1P [6] and G9P [6] in Pakistan [28]. True genotypic diversity of RVA in developing countries like Bangladesh may be greater than we can detect. In future, more studies should be conducted to measure the actual burden of RVA in Bangladesh.
Based on the divergence of VP7 genes of G1 rotaviruses a classification scheme has been developed that specifies G1 rotavirus into various lineages and sublineages. Numerous studies had been undertaken to specify the heterogeneity and dynamics of the evolution of G1 rotaviruses in USA, Italy and Japan [23,38,39]. In this study, we followed the classification scheme that was used in Phan [8] and G2P [4] genotypes in Bangladesh showed close relationship with previously isolated strains in Bangladesh. Our study is an important tool to understand the genetic diversity of circulating rotavirus in children with AGE in Bangladesh. The presence of diverse genotypes with mutation in antigenic regions can cause reinfection in vaccinated population and reduction of vaccine efficiency in Bangladesh. Further, circulation of diverse genotypes is the major reason of high prevalence of rotavirus in Bangladesh [17,40].
Rotavirus VP7 protein has three important antigenic regions-amino acid position 87 to 100 (A), 141-150 (B) and 208-224 (C) which are involved in interaction with antibodies and neutralization of epitopes [41,42]. We detected mutations at antigenic determination regions in G1P [8] and G2P [4] in Bangladesh during 2014-2019. Mutations at antigenic regions may alter rotavirus antigenic properties [41,42]. Prevalence of unusual genotypes, intragenotypic diversity and frequent mutations at antigenic regions of rotavirus A are major explanations of the constant predominance of rotavirus in Bangladesh.
Bangladesh is one of the countries with the highest disease burden of RVA infections in children <5 years of age. The introduction of RVA vaccines has reduced the global disease burden of RVA in many countries [29]. According to recent survey by Mahmud-Al-Rafat et al.
(2018), RVA vaccine coverage in capital Dhaka is good but in small town and suburban areas the RVA vaccine coverage are poor. Overall, RVA vaccine introduction in EPI program in Bangladesh could prevent over 10,000 deaths per year [19]. Besides, RVA studies including large number of samples for sustained time periods and diverse sentinel sites throughout the country will be required for assessing the true RVA disease burden. Moreover, whole genome analysis of human RVA strains and inclusion of animal and environmental samples in future studies is recommended to elucidate the current interspecies transmission and genomic reassortment and recombination events.

Conclusions
This study concludes the high prevalence of G1P [8] and G2P [4] RVA strains along with reduced prevalence of G3 RVA in Bangladesh during 2014-2019. This study will provide comprehensive insights to the researchers and public health authorities to measure the RVA disease burden in the country. The observed large diversity of RVA genotypes along with their yearly and seasonal fluctuations require an in-depth, broad surveillance system in the country. The findings of this study can provide guideline data before introduction of rotavirus vaccine in the national immunization program in Bangladesh. The data in this study can be used to assess the impact of rotavirus vaccine in the future.