Pre COVID-19 molecular epidemiology of respiratory syncytial virus (RSV) among children in Bangladesh.

Epidemiological data of specific respiratory pathogens from the pre-COVID-19 period are essential to determine the effects of the COVID-19 pandemic on other respiratory infections. In this study, we revealed the pre-COVID-19 molecular epidemiology of respiratory syncytial virus (RSV) among children in Bangladesh. We tested 3170 samples collected from 2008 to 2012 for a panel of respiratory viruses; RSV, human metapneumovirus (hMPV), human parainfluenza viruses (hPIV) 1, 2, 3, and adenovirus. Five hundred fifty-five samples (17.5 %) were positive for RSV, including 2.5% having co-infections with other viruses. Genotypic characterization of RSV showed that RSV-A (82%) contributed more acute respiratory infections than RSV-B (18%). Clinical features were similar with RSV-A and RSV-B infections. However, children with RSV-B were more likely to have upper respiratory infections (URI) (10% vs. 29%, p = 0.03). Among RSV-A cases, hospitalization was higher for ON1 cases (25%, ON1 vs. 8%, NA1, p = 0.04), whereas the recovery without a disability was higher among the NA1 cases (56%, ON1 vs. 88%, NA1, p = 0.02). The time to the most recent common ancestor (TMRCA) for RSV in Bangladesh was 1949 for RSV-A and 1944 for RSV-B. This study revealed the genotypic diversity and evolutionary relatedness of RSV strains in Bangladesh and provided pre-COVID molecular epidemiology data to understand better the COVID-19 impact on upcoming RSV epidemiology in Bangladesh.

Introduction

Since the inception of the Coronavirus disease 2019 (COVID-19) pandemic, interventions to control the spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) such as strict lockdowns, limited international movements, social distancing, vaccination, proper hygiene practices, and restricted public gatherings have had indirect effects on the dynamics of other seasonal acute respiratory infections (ARIs), most notably on the RSV infections and influenza [1]. ARIs contribute significant health burden for developing nations such as Bangladesh [2]. During the pre-COVID era, RSV was one of the major causes of ARIs in children under 5 years of age in Bangladesh [3, 4, 5]. About 99% of RSV-associated deaths occur in low-income countries, ten times more than in high-income areas [3, 6]. In 2015, 3.2 million hospitalization and 60 thousand in-hospital deaths occurred globally in children younger than five years of age due to RSV [7]. About 50–70% of infants experience their first RSV infection in the first year of life, and most of them are infected by the age of two [7, 8, 9]. Susceptibility to RSV decreases with age but also increases in the elderly (>64 years of age) [9, 10, 11]. RSV season lasts variably between three to nine months as multiple short epidemics and has no specific pattern [12, 13].

RSV is an enveloped, non-segmented, negative-sense, single-stranded RNA virus that belongs to the genus Orthopneumovirus and family Pneumoviridae [14, 15]. Based on the second hyper-variable region (HVR2) of glycoprotein G, two antigenically distinct subgroups of RSV-A and RSV-B have been identified [16, 17, 18]. The nucleotide sequence variability in the G gene is higher than the other RSV genes and can reveal the genetic diversity of strains [19, 20]. The prevalence of RSV-A and RSV-B can change from year to year [21, 22]. RSV-A has 14 genotypes (ON1-2, GA1-GA7, SAA1, and NA1-4), whereas RSV-B so far has 24 reported genotypes (GB1-GB4, SAB1-SAB4, URU1-2, CB, BA1–BA12, and BAC) [6, 8, 10, 23].

There is limited knowledge on the RSV genotypes and their evolutionary patterns in Bangladesh. Like all RNA viruses, RSV circulates in the community, causing febrile illness, and undergoes genetic variation that can lead to aggressive infections and potential outbreaks [18]. The ongoing COVID-19 pandemic may become seasonal COVID-19 outbreaks with time, consequently impacting the epidemiology of respiratory pathogens worldwide [24]. For example, several southern and northern hemisphere countries have now reported delayed RSV peaks [25]. The objective of this study was to understand the molecular epidemiology and evolutionary relatedness of RSV in relation to clinical presentations among children in Bangladesh during the pre-COVID-19 period. As we hopefully enter a post-pandemic era when RSV vaccines and new therapies are likely to become available, our findings will help understand the impact of COVID-19 on RSV epidemiology in Bangladesh.

Materials and methods

Real-time reverse transcription polymerase chain reaction (RT-PCR) of pre-COVID-19 respiratory samples

A total of 3170 archived nasopharyngeal wash (NPW) specimens, previously tested negative for influenza viruses, were selected for this study (Figure 1 ). The specimens were collected from viral ARI patients of less than five years of age, between 2008 and 2012, under the study “Surveillance for influenza and the viral etiologies of influenza-like febrile illnesses in an urban slum in Dhaka, Bangladesh.” The study was approved by the Research Review Committee (RRC) and Ethical Review Committee (ERC) of International Centre for Diarrhoeal Disease Research (icddr,b) (Protocol no: PR: 2003-030). A field research assistant (FRA) identified patients with reportable illness (fever, cough, fast or difficult breathing, noisy breathing, inability to feed, or lethargy) during their weekly visits to selected households with children <5 years at the Kamalapur surveillance area in Dhaka, Bangladesh. Patients were interviewed using icddr,b institutional review board (IRB) approved questionnaires, and clinical data were documented in case report forms. If fever (≥38 °C axillary) was reported or elevated respiratory rate (RR) documented, then the FRA referred the child to the field clinic for examination by a medical officer (MO). The child’s RR was considered elevated if it was ≥60/minute among children <2 months, ≥ 50/minute among children between 2 and 11 months, and ≥40/minute among children 12–59 months. If the child was identified as suspected viral ARI (fever of ≥38 °C or RR, plus at least one of the following: cough, chest-indrawing, crepitations (inspiration), wheezing (expiratory), or rhonchi) when evaluated by the MO, then he/she was eligible for study enrolment. Parents or legal guardians were informed about the study, and written consent was received before collecting specimens and clinical data. The MO used sterile techniques and precautions to collect NPW samples by flushing and aspirating 3 ml sterile normal saline into one nostril using a 3 ml sterile polypropylene syringe with an attached butterfly assembly without the needle. The samples were collected into sterile containers containing viral transport media and transported to the virology lab at icddr,b in a cool box (4.0 ± 2.0 °C). Aliquots of these archived samples were kept frozen at or below −70 °C, and samples were thawed before use. Viral RNA was extracted from 200 μL NPW samples using InviMag®Virus DNA/RNA Mini Kit/KF 96 (STRATEC Molecular, Germany) kit. Extracted RNA was used as a template for detecting a panel of common respiratory viruses like RSV, HPIV-1, 2, 3, HMPV, and adenoviruses by real-time RT-PCR assays [26, 27]. For further genotypic analysis, every fifth of the solely positive RSV samples (n = 94) was chosen randomly (Figure 1).

Figure 1

Consort flow diagram showing sample disposition throughout the study for the genetic characterization of RSV.

RSV genotyping

RNA from RSV PCR-positive cases was used to prepare complementary DNA (cDNA) using a High Capacity cDNA kit (Applied Biosystems™, Foster City, USA). External and heminested PCRs were performed using the Qiagen Hotstar TaqDNA Polymerase Kit (Qiagen, Hilden, Germany), targeting the HVR2 region of the G gene of RSV according to the procedure described elsewhere [28, 29]. The forward primers used for the first and second rounds of heminested PCR amplification of RSV-A G gene were 5′-GAAGTGTTCAACTTTGTACC-3′ and 5′-TATGCAGCAACAATCCAACC-3′ respectively. For the first and second rounds of RSV-B G gene PCR amplification, the forward primers were 5′-AAGATGATTACCATTTTGAAGT-3′ and 5′-TATGCAGCAACAATCCAACC-3′ respectively. The reverse primer 5′-CAACTCCATTGTTATTTGCC-3′ was used for the two rounds of both RSV-A and RSV-B G gene heminested PCR amplification. After heminested PCR, the products were analyzed through gel electrophoresis on 1.5% agarose gels. The characteristic PCR product size was 458 bp for RSV-A and 464 bp for RSV-B. All RSV positive PCR products were purified using the ExoSAP-IT (Affymetrix, California, USA). The cycle sequencing reaction was performed using the ABI BigDye Terminator v.3.1 cycle sequencing kit (Applied Biosystems, Foster City, USA), and Sanger sequencing was carried out in an ABI 3500 XL genetic analyzer (Applied Biosystems, Foster City, USA) using the consensus forward and reverse primers separately.

Bayesian phylogenetic analyses of RSV-A and RSV-B G gene

To estimate the evolution rate and reconstruct the demographic history of RSV strains, 58 RSV-A and 114 RSV-B nucleotide sequences from GenBank and 58 sequences from this study along with annotated sample collection date and country were included in the Bayesian phylogenetic analyses (Tables 1 and 2 ). Temporal phylogenies of individual RSV genotypes were inferred in BEAST v1.8.324 using a strict clock, the HKY substitution model, and a constant coalescent tree prior. Bayesian Markov Chain Monte Carlo (MCMC) chains were run for 100 million steps with sampling every 10,000 generations. Run convergence was confirmed in Tracer v1.6 after 10% burn-in removal. Maximum clade credibility (MCC) trees were generated using Tree Annotator v1.7. All phylogenetic trees were visualized using FigTree v1.4.2.

Table 1

List of accession numbers of RSV strains isolated from this study.

RSV Strain	Accession Number	RSV Strain	Accession Number
BGD/101812/RSV_A/2008	MW306451	BGD/407306/RSV_A/2012	MW306480
BGD/231567/RSV_A/2008	MW306452	BGD/341072/RSV_A/2012	MW306481
BGD/282293/RSV_A/2008	MW306453	BGD/551097/RSV_A/2012	MW306482
BGD/336905/RSV_A/2008	MW306454	BGD/455575/RSV_A/2012	MW306483
BGD/221420/RSV_A/2010	MW306455	BGD/516783/RSV_A/2012	MW306484
BGD/236003/RSV_A/2010	MW306456	BGD/549336/RSV_A/2012	MW306485
BGD/335327/RSV_A/2010	MW306457	BGD/504344/RSV_A/2012	MW306486
BGD/239928/RSV_A/2009	MW306458	BGD/516784/RSV_A/2012	MW306487
BGD/359479/RSV_A/2010	MW306459	BGD/359990/RSV_A/2012	MW306488
BGD/365974/RSV_A/2010	MW306460	BGD/516789/RSV_A/2012	MW306489
BGD/341471/RSV_A/2010	MW306461	BGD/363288/RSV_A/2012	MW306490
BGD/341504/RSV_A/2010	MW306462	BGD/341589/RSV_A/2012	MW306491
BGD/346416/RSV_A/2010	MW306463	BGD/343143/RSV_A/2012	MW306492
BGD/347260/RSV_A/2010	MW306464	BGD/202110/RSV_B/2011	MW306493
BGD/350396/RSV_A/2010	MW306465	BGD/221460/RSV_B/2011	MW306494
BGD/216099/RSV_A/2010	MW306466	BGD/356609/RSV_B/2011	MW306495
BGD/350734/RSV_A/2010	MW306467	BGD/367144/RSV_B/2011	MW306496
BGD/353808/RSV_A/2010	MW306468	BGD/407333/RSV_B/2011	MW306497
BGD/354688/RSV_A/2010	MW306469	BGD/407698/RSV_B/2011	MW306498
BDG/343226/RSV_A/2009	MW306470	BGD/516103/RSV_B/2011	MW306499
BDG/357189/RSV_A/2010	MW306471	BGD/200050/RSV_B/2008	MW306500
BGD/350670/RSV_A/2010	MW306472	BGD/212919/RSV_B/2008	MW306501
BGD/353848/RSV_A/2010	MW306473	BGD/333335/RSV_B/2008	MW306502
BGD/479069/RSV_A/2011	MW306474	BGD/358743/RSV_B/2011	MW306503
BGD/363369/RSV_A/2009	MW306475	BGD/359943/RSV_B/2012	MW306504
BGD/345145/RSV_A/2012	MW306476	BGD/360300/RSV_B/2012	MW306505
BGD/355224/RSV_A/2012	MW306477	BGD/364511/RSV_B/2011	MW306506
BGD/518783/RSV_A/2012	MW306478	BGD/364699/RSV_B/2012	MW306507
BGD/357322/RSV_A/2012	MW306479	BGD/516950/RSV_B/2011	MW306508

Table 2

List of RSV-A and RSV-B GenBank reference sequences.

RSV-A RefSeq	Country	Collection Date	RSV-A RefSeq	Country	Collection Date
KF246605	INDIA	2010	JX627337	CUBA	2012
JN257699	CANADA	2010	JX627338	USA	2012
JX513282	BRAZIL	2009	JX627339	SPAIN	2012
JF920052	USA	2010	JX627340	INDIA	2011
KT781346	CHINA	2014	JX627341	INDIA	2012
KJ627304	PERU	2009	JX627342	PHILIPPINES	2012
KC476706	SOUTH AFRICA	2010	JX627343	PARAGUAY	2013
JX015496	NETHERLANDS	2008	JX627344	INDIA	2012
KP792352	NETHERLANDS	2007	JX627345	KENYA	2012
JX015481	NETHERLANDS	2009	JX627346	SPAIN	2012
JX015482	NETHERLANDS	2006	JX627347	CUBA	2010
KF246618	INDIA	2010	JX627348	SPAIN	2015
KF826821	USA	2007	JX627349	SOUTH AFRICA	2012
JX015492	NETHERLANDS	2007	JX627350	CHINA	2014
KF826838	ARGENTINA	2006	JX627351	CANADA	2010
KF530261	GERMANY	2008	AF065256	USA	1993
JQ901452	NETHERLANDS	2001	Z33424	SPAIN	1988
JX131637	SAUDI ARABIA	2008	JF920062	USA	1998
KC476743	SOUTH AFRICA	2009	JX069800	USA	1997
KF826848	AUSTRALIA	2007	Z33455	SPAIN	1992
KP317944	KENYA	2006	AF06522	USA	1993
KP317948	KENYA	2005	HQ699263	KOREA	2009
JX069799	USA	2001	JX015485	NETHERLANDS	2005
KJ627260	PERU	2008	KP792353	NETHERLANDS	2007
KJ627327	PERU	2007	KC476656	SOUTH AFRICA	2010
AB754590	JAPAN	2012	AF065254	USA	1993
KM434001	CHINA	2013	Z33431	SPAIN	1991
JX627336	KOREA	2011

Statistical analysis

The demographic and clinical characteristics along with patient outcomes were summarized using frequency and percentage values. To compare the demographic and clinical characteristics of children with different RSV subtypes (RSV-A vs RSV-B and RSV-A ON1 vs NA1 genotype), student’s t-tests were used to compare mean differences. The difference in proportions were assessed using the chi-squared χ2 or Fisher exact test as appropriate. The multivariable Poisson regression adjusted for age and sex analysis was performed to estimate the relative risk and risk difference for those clinical characteristics that showed a significant difference between children with RSV-A ON1 and NA1 infections. In the Poisson regression analysis, RSV-A NA1 ARI group (n = 26) was used as the reference group. For the multivariate regression, variables including age and sex were adjusted for estimating the adjusted Relative Risk (Adj RR) and 95% Confidence Interval (CI). We considered p < 0.05 as statistically significant. Stata 15 (Stata Corp. 2013. Stata Statistical Software: Release 13. College Station, TX: Stata Corp LP.) was used for all analysis.

Results

Out of 3170 specimens, 555 (17.5%) were positive for RSV by real-time RT-PCR, of which 470 (14.8%) were positive only for RSV, and 85 (2.7%) had co-infection of RSV and other viral pathogens included in the testing panel. Adenovirus (13%) was the second most frequently detected after RSV, followed by hPIV (9%) and hMPV (6%). However, 5.7% of the cases (n = 177) were due to co-infections by respiratory viruses, other than RSV, like HPIV-1, 2, 3, HMPV or adenoviruses. The year-wise prevalence of RSV ranged widely from 2.3% (2009) to 24.5% (2010).

Seasonality of RSV

The RSV infections were primarily detected in hot and rainy seasons, which extended from mid-June to October in 2008, 2011, and 2012 except 2009–2010, when the season initiated at the end of 2009 and lasted the first three months of 2010 (Figure 2 ). RSV contributed more viral ARI cases in July–August (36% and 58.4%) in 2008, November–December (21% and 17%) in 2009, January–February (68% and 64%) in 2010, August–September (38% and 52%) in 2011 and September–October (69% and 61%) in 2012. However, there were some pauses in monthly sample collection due to the lack of patients identified with viral ARI symptoms, which may have confounded this seasonality data. The 2009 global swine flu pandemic caused by Influenza A(H1N1)pdm09 virus might have also affected RSV seasonality during the study period, especially in 2009–2010.

Figure 2

Seasonal distribution of RSV detected between 2008 and 2012 in this study.

Prevalence of RSV-A and RSV-B subgroups

Every fifth (20%, n = 94) of RSV PCR-positive samples were further characterized into the RSV-A and RSV-B subgroups. During the study period, both RSV-A (n = 77, 81.9%) and RSV-B (n = 17, 18.1%) were detected by PCR. All (n = 94) samples were subjected to Sanger sequencing to obtain RSV G gene fragment sequences. Due to the low nucleic acid content of some PCR positive RSV samples, band visualization in gel electrophoresis and subsequent retrieval of conclusive RSV nucleotide sequence data was difficult. Therefore, we could only obtain 42 and 16 RSV G gene fragment sequences from RSV-A and RSV-B subgroups, respectively. These sequences were submitted to GenBank (accession no. MW306451-MW306508) and subjected to phylogenetic analysis.

Genetic comparison of RSV-A and RSV-B strains

All RSV genotypes were identified through nucleotide similarity searches performed using the National Center for Biotechnology Information (NCBI, National Institutes of Health, Bethesda, MD) BLAST (Basic Local Alignment Search Tool) server on GenBank database. The genotypes were further confirmed maximum-likelihood (ML) phylogenetic analysis in MEGA 7 using 1000 bootstrap replicates. The nucleotide sequence analysis of the G gene fragment of 42 RSV-A strains showed that 16 (MW306477-MW306492) sequences were ON1 genotypes and 26 (MW306451-MW306476) were NA1. The sequences of the ON1 genotype had the characteristic 72 nucleotides (24 amino acids at 260 to 307 of M11486 prototype RSV strain) duplication in their G gene fragment (Figure 3 ). The RSV strains of the NA1 genotype were most prevalent till 2011 in Bangladesh and gradually decreased after that (Figure 4 ). The strains of the ON1 genotype appeared in the year 2012. The nucleotide and deduced amino acid sequence variations between the ON1 study strains and prototype ON1 strain were up to 2.5% and 4.3%, respectively. The nucleotide and amino acid variations between the NA1 study strains and prototype NA1 strain were up to 6% and 15%, respectively. The nucleotide sequence analysis of RSV-B sequences demonstrated the circulating strains belonged to the GB13 genotype containing the characteristic 60 nucleotides (20 amino acids) duplication in the G gene HVR2 (Figure 5 ). The pairwise distances between strains of RSV-B ranged from 0.1 to 7.3% at the nucleotide level and 0.4 to 10.8% at the virtually transcribed amino acid sequences.

Figure 3

Deduced amino acid alignment and changes in the second hypervariable region of G protein of the NA1 and ON1 genotype. The figure includes the alignment of Bangladeshi RSV-A strains with the prototype strain (M11486), NA1 reference strain (NC_001803), ON1 reference strain 1 (AEQ98758), and ON1 reference strain 2 (AHG54506). The amino acids sequence alignment corresponds to 212–320 amino acids of the prototype strain. Dashes indicate identical residues. Rectangular boxes indicate the two copies of the duplicated 24-amino-acid region in group ON1 strains. Potential N-glycosylation sites are indicated by grey shading.

Figure 4

The seasonal distribution of RSV-A (NA1 and ON1) and RSV-B (GB13) genotypes detected in this study.

Figure 5

Deduced amino acid alignment and changes in the second hypervariable region of G protein of the GB13 genotype BA4 lineage. The figure includes alignment of Bangladeshi RSV-B strains with the prototype strain from Argentina (AY333364). The amino acids sequence alignment corresponds to 207–312 amino acids of the prototype strain. Dashes indicate identical residues. Rectangular boxes indicate the two copies of the duplicated 20-amino-acid region in group BA4 strains. Potential N-glycosylation sites are indicated by grey shading.

Demographic and clinical characteristics of RSV subtypes

The predominant clinical presentation of patients from both RSV-A and RSV-B groups was fever. Clinical features and duration of illness did not vary between children with RSV-A and RSV-B (Median days 2.0 vs. 3, P = 0.55) infections. The majority of RSV-A and RSV-B patients were prescribed with amoxicillin antibiotics (81% vs. 82%, p = 0.95). Compared to children with RSV-A, children with RSV-B were more likely to have a diagnosis of upper respiratory infection (URI) (10% vs. 29%, p = 0.03). The majority of children with both RSV-A and RSV-B recovered (77% vs. 76%, p = 0.98) from their illness (Table 3 ). Among the RSV-A patients, hospitalization was higher among the ON1 ARI patients (25%, ON1 vs. 8%, NA1, p = 0.04), whereas the recovery without a disability was higher among the NA1 ARI patients (88%, NA1 vs. 56%, ON1 p = 0.02) (Table 4 ). The ON1 ARI patients had a higher risk of difficulty in breathing (Adj RR: 1.74 (95% CI: 1.1–2.5)), tachycardia (Adj risk difference: 8.7 (95% CI: 1.1–16.4)) and tachypnea (Adj risk difference: 4.5 (95% CI: 0.0–9.1)) after adjusting demographic confounders. Though they were significant in univariate analysis, other clinical features such as fever, hospitalization, and treatment outcome did not differ in multivariate analysis (Tables 4 and 5 ).

Table 3

Demographic and clinical characteristics of children with RSV-A and RSV-B, Dhaka, 2008–2012.

	RSV-A (n = 77)	RSV-B (n = 17)	p value
Demographic information
Sex (male) (%)	44 (57.14)	9 (52.94)	0.752
Median age (in months	11 (IQR:8–21)	13 (IQR:8–24)	0.764
Age group
0–12 month	42 (54.55)	9 (52.94)	0.667
13–24 month	19 (24.68)	4 (23.53)
25–36 month	10 (12.99)	4 (23.53)
37–48 month	3 (3.90)	-
49–60 month	3 (3.90)	-
Clinical information
Chief complaint			0.898
Cough	10 (12.99)	2 (11.76)	0.890
Difficult breathing	28 (36.36)	5 (29.41)	0.586
Fast breathing	6 (7.79)	1 (5.88)	0.786
Fever	33 (42.86)	9 (52.94)	0.449
Duration of chief complaint	2 (IQR:2–4)	3 (IQR:2–3)	0.553
Fever (%)	69 (89.61)	17 (100)	0.164
Headache (%)	56 (72.73)	14 (82.35)	0.41
Eye pain (%)	57 (74.03)	14 (82.35)	0.47
Vomiting (%)	76 (98.70)	16 (94.12)	0.236
Runny nose (%)	66 (85.71)	14 (82.35)	0.724
Cough (%)	76 (98.70)	17 (100)	0.636
Difficult breathing (%)	49 (63.64)	7 (41.18)	0.087
Fast breathing (%)	51 (66.23)	8 (47.06)	0.138
Ear pain (%)	1 (1.30)	0	0.636
Ear discharge (%)	1 (1.30)	0	0.636
Temperature (°C) (mean (95% CI))	37.95 (95% CI: 37.73–38.17)	38.43 (95% CI: 38.11–38.75)	0.052
Pulse (per minute) (mean (95% CI))	136.68 (95% CI: 133.59–139.78)	142.41 (95% CI: 134.32–150.49)	0.138
Respiratory (%) (mean (95% CI))	52.25 (95% CI: 50.48–54.03)	52.70 (95% CI: 48.38–57.02)	0.836
Spo 2 (mean (95% CI))	96.90 (95% CI: 96.60–97.21)	97.11 (95% CI: 96.26–97.96)	0.589
Ears			0.798
Normal	75 (97.40)	17 (100)
Pain	1 (1.30)	-
Other abnormality	1 (1.30)	-
Nose			0.904
Normal	35 (45.45)	8 (47.06)	0.904
Clearly discharge	42 (54.55)	9 (52.94)	0.904
Pharynx	77 (100)	17 (100)
Lymph nodes	77 (100)	17 (100)
Respiratory auscultation			0.309
Clear	14 (18.18)	6 (35.29)	0.118
Crepitation	4 (5.19)	0	0.337
Crepitation with wheezes	59 (76.62)	11 (64.71)	0.308
Treatment and outcome
Hospitalized			0.317
Yes	10 (12.99)	0
No	64 (83.12)	16 (94.12)
Do not know	3 (3.90)	1 (5.88)
Preliminary clinical diagnosis			0.208
Pneumonia	4 (5.19)	0	0.337
Pneumonia with wheeze	51 (66.23)	10 (58.82)	0.562
Severe pneumonia with wheeze	8 (10.39)	1 (5.88)	0.567
Typhoid fever	4 (5.19)	4 (23.53)	0.014
Upper Respiratory Illness	10 (12.99)	2 (11.76)	0.89
Final diagnosis			0.56
Bronchiolitis	1 (1.30)	0	0.636
Pneumonia	3 (3.90)	0	0.407
Pneumonia with wheeze	50 (64.94)	10 (58.82)	0.634
Severe pneumonia with wheeze	13 (16.88)	2 (11.76)	0.601
Sinusitis	1 (1.30)	0	0.636
Typhoid fever	1 (1.30)	0	0.636
Upper Respiratory Illness	8 (10.39)	5 (29.41)	0.039
Outcome of illness			0.862
Recovered	59 (76.62)	13 (76.47)	0.989
Recovered with disability	15 (19.48)	4 (23.53)	0.706
Unknown	3 (3.90)	0	0.407
Duration of illness (mean (95% CI))	8.01 (95% CI: 6.71–9.31)	7.64 (95% CI: 4.85–10.45)	0.813

Table 4

Demographic and clinical characteristics of children with RSV-A ON1 and NA1 genotype, Dhaka, 2008–2012.

	RSV-A ON1 (n = 16)	RSV-A NA1 (n = 26)	p value
Demographic information
Sex (male) (%)	6 (37.50)	18 (69.23)	0.043
Median age (in month)	9.36 (IQR:8.03–21.69)	17.19 (IQR:12.16–26.13)	0.15
Age group
0–12 month	12 (75)	8 (30.77)	0.026
13–24 month	0	10 (38.46)
25–36 month	2 (12.50)	6 (23.08)
37–48 month	1 (6.25)	1 (3.85)
49–60 month	1 (6.25)	1 (3.85)
Clinical information
Chief complaint			0.026
Cough	0	3 (11.54)	0.159
Difficult breathing	10 (62.50)	6 (23.08)	0.012
Fast breathing	3 (18.75)	3 (11.54)	0.524
Fever	3 (18.75)	14 (53.85)	0.026
Duration of chief complaint	2 (IQR:1.5–3.5)	2 (IQR:2–4)	0.239
Fever (%)	11 (68.75)	24 (92.31)	0.05
Headache (%)	16 (100)	14 (53.85)	0.001
Eye pain (%)	16 (100)	14 (53.85)	0.001
Vomiting (%)	16 (100)	26 (100)	-
Runny nose (%)	15 (93.75)	24 (92.31)	0.862
Cough (%)	16 (100)	26 (100)	-
Difficult breathing (%)	15 (93.75)	14 (53.85)	0.006
Fast breathing (%)	15 (93.75)	14 (53.85)	0.006
Temperature (°C) (mean (95% CI)	37.5 (95% CI: 37.10–37.89)	37.86 (95% CI: 37.45–38.26)	0.233
Pulse (per minute) (mean (95% CI))	140.75 (95% CI: 134.69–146.80)	131.53 (95% CI: 127.09–135.97)	0.015
Respiratory (%) (mean (95% CI))	53.75 (95% CI: 50.89–56.60)	49.38 (95% CI: 46.37–52.39)	0.054
Spo 2 (mean (95% CI))	96.93 (95% CI: 95.98–97.88)	96.92 (95% CI: 96.34–97.50)	0.978
Eyes	16 (100)	26 (100)	-
Normal	16 (100)	25 (96.15)
Pain	0	1 (3.85)
Nose			0.85
Normal	6 (37.50)	9 (34.62)	0.852
Clearly discharge	10 (62.50)	17 (65.38)	0.852
Pharynx	16 (100)	26 (100)	-
Lymph nodes	16 (100)	26 (100)
Respiratory auscultation			0.033
Clear	0	7 (26.92)	0.023
Crepitation	0	1 (3.85)	0.427
Crepitation with wheezes	16 (100)	18 (69.23)	0.014
Treatment and out come
Hospitalized			0.043
Yes	4 (25)	2 (7.69)
No	10 (62.50)	24 (92.31)
Do not know	2 (12.50)	0
Preliminary clinical diagnosis			0.08
Pneumonia	0	1 (3.85)	0.427
Pneumonia with wheeze	13 (81.25)	17 (65.38)	0.274
Severe pneumonia with wheeze	3 (18.75)	1 (3.85)	0.12
Typhoid fever	0	2 (7.69)	0.256
Upper Respiratory Illness	0	5 (19.23)	0.062
Final diagnosis			0.071
Bronchiolitis	0	1 (3.85)	0.427
Pneumonia	0	1 (3.85)	0.427
Pneumonia with wheeze	11 (68.75)	19 (73.08)	0.766
Severe pneumonia with wheeze	5 (31.25)	2 (11.76)	0.127
Sinusitis	1 (1.30)	1 (3.85)	0.633
Upper Respiratory Illness	0	3 (11.54)	0.159
Outcome of illness			0.037
Recovered	9 (56.25)	23 (88.46)	0.02
Recovered with disability	5 (31.25)	3 (11.54)	0.122
Unknown	2 (12.50)	0	0.075
Duration of illness (mean (95% CI))	7.87 (95% CI: 4.68–11.06)	7.11 (95% CI: 5.27–8.96)	0.656

Table 5

Multivariate Poisson regression analysis of clinical severity comparison between RSV-A ON1 and NA1 genotype ARI cases.

RSV subgroup A positive case (N = 42)
	RSV-A ON1 (n = 16)	∗RSV-A NA1 (n = 26)(reference group)	P value	Unadjusted RR	95% CI	∗Adjusted RR	95% CI
Clinical information
Difficult breathing (%)	15 (93.75)	14 (53.85)	0.012	1.74	1.2–2.5	1.74	1.1–2.5
Fever (%)	11 (68.75)	24 (92.31)	0.05	0.74	0.36–1.5	0.72	0.34–1.5
Pulse (per minute) (mean (95% CI)) (risk difference)	140.75 (95% CI: 134.69–146.80)	131.53 (95% CI: 127.09–135.97)	0.015	9.2	2.0–16.3	8.7	1.1–16.4
Respiratory (%) (mean (95% CI)) (risk difference)	53.75 (95% CI: 50.89–56.60)	49.38 (95% CI: 46.37–52.39)	0.054	4.3	0.04–8.6	4.5	0.0–9.1
Treatment and outcome
Hospitalization	4 (25)	2 (7.69)	0.043	1.02	0.66–1.5	1.05	0.67–1.6
Outcome of illness (recovered)	9 (56.25)	23 (88.46)	0.02	0.72	0.33–1.5	0.72	0.32–1.62

In the Poisson regression analysis, RSV-A NA1 ARI group (n = 26) was used as the reference group. In the multivariate regression, variables including age and sex were adjusted for estimating adjusted Relative Risk (Adj RR) and 95% Confidence Interval (CI).

Phylodynamics and evolutionary relationships of RSV strains

The effective sample size (ESS) value after the BEAST run was 457 and 123 for RSV-A and B respectively (Figure 6 ). The time-scaled MCC tree illustrated the time to the most recent common ancestor (TMRCA) of the Bangladeshi RSV strains, which was estimated to be 1949 considering 95% highest probability density (HPD): 1943–1954 for RSV-A (Figure 7 ). The onset of genotype GA2 was around 1980 [95% HPD 1975–1985]. After that, GA2 became the parent genotype for the evolution of NA1 and ON1 in 1996 [95% HPD 1993–1998] and 2007 [95% HPD 2004–2010], respectively. For RSV-B, the TMRCA was estimated to be ∼1944 [95% HPD: 1934–1954] (Figure 8 ). The molecular evolutionary rate was 2.07 × 10−3 substitutions/site/year [95% HPD: 1.65–2.51 × 10−3] for RSV-A strains and 2.34 × 10−3 substitutions/site/year [95% HPD: 1.92–2.76 × 10−3] for RSV-B strains.

Figure 6

Summary statistics and posterior distribution frequency graph for BEAST analysis of A RSV-A and B RSV-B using HKU substitution model with 100 million iterations and sampling every 10,000 generations.

Figure 7

Time-scaled Bayesian maximum clade credibility (MCC) trees for HRSV-A G gene sequences.

Figure 8

Time-scaled Bayesian maximum clade credibility (MCC) trees for HRSV-B G gene sequences.

Amino acid sequence and N-linked glycosylation site analysis

Multiple sequence alignment of the deduced amino acid sequence of the HVR2 of Bangladeshi RSV-A strains was performed with the prototype RSV-A strain A2 (M11486), NC_001803, and reference strains from India (AHG54430, AHG54512, AHG54506) and Canada (AEQ98758). The predicted amino acid sequences of RSV-A strains obtained in this study corresponded to residue position 198–321 (based on RSV-A strain A2). Amino acid alignment revealed that a 24 residue duplication was found in the ON1 sequence. Several infrequent amino acid substitutions were found in the characteristic 24 residue duplication of Bangladeshi ON1 strains (Figure 3). Compared to RSV-A reference and prototype strains, the Bangladeshi RSV NA1 strains had N→D and P→L substitutions at 237 and 274 positions, respectively. Such substitutions were not observed among the ON1 strains.

Nevertheless, several infrequent amino acid substitutions in both the NA1 and ON1 genotypes revealed the strain diversity of RSV-A in Bangladesh. Four potential N-linked glycosylation sites at positions 237, 251, 273, and 318 were identified in the partial amino acid sequences of RSV-A G protein. N-linked glycosylation was most conserved at 318 position in all ON1 and NA1 strains. Most of the NA1 strains lost the glycosylation site at 237 due to the N→D substitution. Among the strains of the ON1 genotype, the N-linked glycosylation site was conserved at positions 237 and 251 but was altered at 273 through the substitution of N→Y (Figure 3).

All Bangladeshi RSV-B strains were of GB13 genotype and had the T229I and S247P substitutions that specified a BA4 lineage. Most of them (14/16, 87.5%) also contained the BA4 H287Y substitution and the BA2 L219P substitution (Figure 5). Two potential N-linked glycosylation sites, 230 and 296 (based on BA4128/99B numbering) of RSV-B G protein, remained conserved in Bangladeshi GB13 genotype strains (Figure 5). The Bangladeshi GB13 strains were genetically diverse compared to the reference strain, but some consistent amino acid substitutions were observed.

Discussion

We tested respiratory samples collected on average ten years before the first case of COVID-19 was identified in Bangladesh for retrospective analysis. Out of 3170 samples, 555 (17.5 %) were positive for RSV, including 85 (2.5%) samples having co-infections with other respiratory viruses like hMPV, hPIV 1, 2, 3, and adenovirus. Genetically diverse strains of RSV were circulating in Bangladesh during 2008–2012, and molecular analysis of HVR2 of the G protein revealed that the study strains belonged to RSV-A (NA1 and ON1) or RSV-B (GB13) genotypes.

According to our clinical findings, children with RSV-A infections were hospitalized more than children with RSV-B infections, though the clinical manifestations were similar irrespective of the genotypes (Table 3). RSV-A ON1 strains are known to cause URIs more than other RSV-A genotypes like NA1 and GA2 [32,33]. No statistically significant differences were found in the clinical features between ON1 and NA1 subtypes (Table 4). However, the attribution of NA1 in URIs was relatively higher compared to ON1 strains. Our findings suggest that hospitalization due to RSV infections may be related to environmental or host factors.

The phylogenetic analysis of G gene sequences revealed that Bangladeshi RSV-A strains clustered into two genotypes: NA1 and ON1 (Figure 7). The NA1 genotype was prevalent in the consecutive years of 2008, 2009, and 2010. The NA1 strains generally possess either the 237 or 273 glycosylation site [10]. This phenomenon was also observed in this study, where the glycosylation at 273 was lost in NA1 strains (Figure 3). In 2011, the ON1 genotype was present and had the characteristic 72 nucleotide insertion (24 amino acid) [17, 30]. Our analysis revealed a pairwise distance of 0.025 at the nucleotide level between study strains and the ON1 prototype strain, which falls within the recommended range (p-distance < 0.049) for clade designation [31]. RSV-B GB13 strains contained the characteristic H287Y substitution specific for BA4 lineage and the BA2 associated L219P substitution (Figure 5). These results match previous reports from various countries [32, 33, 34]. The 72 nucleotide insertion for RSV-A and 60 nucleotide insertion for RSV-B in the G gene may contribute to RSV’s ability to re-infect individuals [35].

Based on our BEAST analysis, the TMRCA of RSV-A and RSV-B was estimated to be 1949 and 1944, respectively. For NA1 and ON1 genotype, TMRCA was estimated to be 1996 and 2007, respectively. Previous reports estimated that the phylogenetic branching time for RSV-A was in 1947 [95% HPD 1930–1956], for RSV-B in 1953 [95% HPD 1938–1962], for NA1 genotype in 1998 [95% HPD 1996–1999] and for ON1 genotype in 2005 [95% HPD 2003–2006] [36,37]. NA1 is a known ancestor of the ON1 genotype and was first detected in Japan in 2004. The first ON1 strain was reported in Canada in 2010–2011 [9,38]. In this study, the ON1 strains were also found in samples collected in 2012, and our data indicated that the ON1 genotype evolved for five years before its first reported case. During 2008–2012, RSV-A and RSV- B were co-circulating. The molecular evolutionary rate of RSV-A was previously reported to be faster than that of RSV-B strains [37, 39]. According to our analysis, the estimated mean evolutionary rate of Bangladeshi RSV-B strains (2.34 × 10−3 substitutions/site/year) was mathematically higher than that of RSV-A strains (2.07 × 10−3 substitutions/site/year), but they did not differ significantly. However, the values were close to previously reported estimates of RSV-A (1.83 × 10−3 to 4.68 × 10−3 substitutions/site/year) and RSV-B (1.95 × 10−3 to 5.89 × 10−3 substitutions/site/year) evolutionary rates [32]. ESS values above 200 indicate confident BEAST analysis results; this was achieved for the RSV-A run only [40]. Due to the limited number of evaluable RSV-B sequences in this study, the ESS value was low. Though larger ESS values are better, Tracer flags up ESS <100. As we conducted an opportunistic research with archived samples for retrospective analysis the sample size could not be increased in the present day. For future research, the low ESS value after RSV-B BEAST run may be improved by using a different model or increasing iterations.

This study was conducted using selected archived samples that were collected during the pre-COVID-19 period. The study was limited to the Kamalapur community at Dhaka city and there was a lack of sample collection during off-season months as no suspected viral ARI patients were identified. It might have confounded the conclusion on RSV seasonality which considered RSV positive samples only. The low nucleic acid content in some of the PCR positive samples made the detection of PCR amplified bands in gel electrophoresis and subsequent sequence analysis difficult. As a result, the small number of RSV sequences evaluated might not provide a true representation burden of infections in Bangladesh from 2008 to 2012. Due to the unavailability of RSV culture facilities and lack of additional funding to pursue viral culture services from contract research organizations, the next-generation sequencing and the complete genomic profiling of RSV strains were not possible.

Many researchers believe that COVID-19 is likely to become a seasonal disease as populations achieve herd immunity [24]. This may impact and revise the seasonality and epidemiology of other respiratory pathogens like RSV [1]. Emerging SARS-CoV-2 variants are reportedly more transmissible and infect humans from a broader age group [41, 42]. Children tend to experience milder COVID-19 symptoms than adults, and the death rate in children is low, so public health interventions against COVID-19 like vaccination are less urgent for them [43, 44]. However, children are likely to spread SARS-CoV-2, and this may exert selective pressure in children for COVID-19 infections and transmissions. This study was designed to generate pre-COVID molecular epidemiology data for RSV infections among children in Bangladesh. Our findings will provide baseline data of RSV infections among children in Bangladesh to assess how the COVID-19 pandemic affected RSV epidemiology, seasonality, vaccine design, and implementation. Further comprehensive research involving larger patient groups in both hospital and community settings will improve the understanding of the evolutionary trajectory of emerging genotypes of RSV and its molecular epidemiology from the intra and post-pandemic timeframes.

Institutional review board statement

The study protocol was approved by icddr,b’s institutional review board (IRB) (protocol no: PR: 2003-030).

Informed consent statement

The parents or legal guardians of all study participants gave written informed consent before specimen and data collection.

Declarations

Author contribution statement

Mohammad Enayet Hossain: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper.

Mohammed Ziaur Rahman: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Md. Muzahidul Islam, Ananya Ferdous Hoque, Mariya Kibtiya Sumiya: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Mst. Noorjahan Begum, Mohammad Mamun Alam, K. M. Main Uddin: Analyzed and interpreted the data; Wrote the paper.

Md. Zakiul Hassan: Analyzed and interpreted the data.

W. Abdullah Brooks: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data.

Mustafizur Rahman, Doli Rani Goswami: Contributed reagents, materials, analysis tools or data.

Funding statement

The work was funded by the (CDC) (PI: Abdullah Brooks, GR00720, PR: 2003-030, CoAg No. 5U01CI000628).

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.