Genetic variability human respiratory syncytial virus subgroups A and B in Turkey during six successive epidemic seasons, 2009-2015.

Human respiratory syncytial virus (HRSV) is most important viral respiratory pathogen of acute lower respiratory tract infections in infants and young children worldwide. The circulating pattern and genetic characteristics in the HRSV attachment glycoprotein gene were investigated in Turkey during six consecutive seasons from 2009 to 2015. HRSVA was dominant in the all epidemic seasons except 2011-2012 season. Partial sequences of the HVR2 region of the G gene of 479 HRSVA and 135 HRSVB were obtained. Most Turkish strains belonged to NA1, ON1, and BA9, which were the predominant genotypes circulating worldwide. Although three novel genotypes, TR-A, TR-BA1, and TR-BA2, were identified, they were not predominant. Clinical data were available for 69 HRSV-positive patients who were monitored due to acute lower respiratory tract illness. There were no significant differences in the clinical diagnosis, hospitalization rates, laboratory findings and treatment observed between the HRSVA and HRSVB groups, and co-infections in this study. The major population afflicted by HRSV infections included infants and children between 13 and 24 months of age. We detected that the CB1, GB5, and THB strains clustered in the same branch with a bootstrap value of 100%. CB-B and BA12 strains clustered in the same branch with a bootstrap value of 65%. The BA11 genotype was clustered in the BA9 genotype in our study. The present study may contribute on the molecular epidemiology of HRSV in Turkey and provide data for HRSV strains circulating in local communities and other regions worldwide.

INTRODUCTION

HRSV is the leading cause of lower respiratory tract infections (LRTI), such as bronchiolitis and pneumonia, in infants and children under 5 years of age. It causes morbidity and mortality in high‐risk populations (premature infants, children with heart and lung disease, immunocompromised individuals, etc.).1, 2 HRSV is also a significant cause of morbidity and mortality in adults. It has an impact approaching that of non‐pandemic influenza virus. Reinfection with HRSV can occur throughout life, and annual epidemics occur during winter months in temperate regions.3

HRSV (Family: Paramyxoviridae, Ordo: Mononegavirales) is an enveloped virus with a single‐stranded, negative sense, non‐segmented RNA genome that is 15.2 kb in length. The viral genome consists of 10 genes, coding for 11 proteins in the following order: NS1, NS2, N, P, M, SH, G, F, M2‐1, M2‐2, and L. Two groups, HRSVA and HRSVB, have been described based on reactions with monoclonal antibodies against the G and F glycoproteins.4, 5 The G protein is a type‐II integral membrane protein consisting of 289‐299 amino acids that facilitates virus attachment to a susceptible cell host. The ectodomain of the G protein consists of two hypervariable regions, HVR1 and HVR2, that are split by a highly conserved 13‐amino acid length (aa 155–206) domain. The addition of O‐ and N‐oligosaccharides to G glycoprotein changes its antigenic characteristics.1, 4, 6, 7, 8 HVR2 region, located at the C‐terminal end of protein, commonly used for molecular analyses. It has highest degree of diversity between subgroups and within each subgroup.4, 5, 6 Until now, considering both antigenic and genetic analysis, 14 genotypes of HRSVA (GA1‐GA7, SAA1, CB‐A, NA1‐NA4, and ON1)1, 8, 9, 10, 11, 12, 13 and 28 genotypes of HRSVB (BA1‐BA13, GB1‐GB5, SAB1‐SAB4, URU1‐2, CB1, BA‐C, CB‐B, and THB) have been identified.11, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24 The genetic diversity between HRSV genotypes and occurrence of new genotypes may change its pathogenicity and facilitate evasion from existing host immune response.1, 10, 13

Recently, a novel HRSVA genotype, ON1, with a 72‐nucleotide duplication within HVR2 was reported in Canada and disseminated globally.13, 21, 23, 25, 26, 27 Surveillance on HRSV infection outbreaks and molecular analysis of the virus are needed to investigate whether the novel genotypes spread worldwide and can cause larger outbreaks. No data regarding the molecular epidemiology of HRSV have been generated in Turkey thus far. Therefore, the aim of this study is to investigate the circulation patterns of the virus and to conduct phylogenetic and molecular analysis of different genotypes of HRSVA and B strains from 2009 to 2015.

MATERIALS AND METHODS

Specimen collection and viral RNA isolation

Sixty‐nine nasopharyngeal aspirates (NPAs) were collected from infants or hospitalized children with acute lower respiratory tract infection (ALRTI) from the Paediatric Allergy and Immunology Clinic, Department of Children's Health Paediatric Haematology Oncology Training and Research Hospital, Ankara, from January 2009 to December 2013. Randomly collected samples from other Paediatric Clinics and from hospitalized adults were also included in this study. The diagnosis of ALRTI was based on the presence of symptoms/signs of acute respiratory illness (cough, fever, nasal discharge, chest pain, sputum production, increased respiratory rate for age, retractions, dyspnea, ronchi or crepitations on auscultation) or radiologic evidence of ALRTI at the time of admission. The clinical features of patients and laboratory findings were evaluated. NPA samples were placed in 2‐3 mL of viral transport medium (minimum essential medium Eagle [MEME], Sigma, Taufkirchen, Germany), vortexed and centrifuged. The supernatants were frozen and stored at −80°C until further testing.

Viral genomic RNA (gRNA) was extracted from 400 μL of NPA supernatant using Qiagen EZ1 Virus Mini Kit v2.0 and EZ1 Advanced extraction robot (Qiagen, Hilden, Germany) according to the manufacturer's instructions with an internal control.

Molecular diagnosis

Each sample was tested for the presence of human respiratory syncytial virus (HRSV) using Fast Track Diagnostics/Respiratory Pathogens 21 commercial kit (Esch‐sur‐Alzette, Luxembourg) (influenza virus A, H1N1pdm, influenza virus B, rhinovirus, respiratory syncytial virus A/B, parainfluenza virus 1/2/3/4, coronavirus OC43/229E/NL63/HKU1, parechovirus, enterovirus, adenovirus, human bocavirus, human metapneumovirus, and IC) using an ABI 7500 real‐time PCR platform (Applied Biosystems, Foster City, CA). RT‐PCR was performed using the manufacturer's protocol (Fast Track Diagnostics/Respiratory Pathogens 21).

HRSV‐positive samples (ct values were between 16 and 32) were then subjected to RT semi‐nested PCR to partially amplify the G protein coding region of gRNA for subsequent phylogenetic analysis.

DNA sequencing and phylogenetic analysis

For genetic analysis of the HRSVs, the carboxy terminus of the G protein coding region was amplified by using the Qiagen One Step RT‐PCR Kit (Qiagen) according to the manufacturer's instructions with ABG490 and F164 primers.28 The expected amplification product sizes for the A, B, and BA (Buenos Aires [BA]) strain that contained sequences of HRSV with a 60‐nucleotide duplication in the attachment [G] protein gene) HRS viruses were 607 bp, 610 bp, and 670 bp, respectively.

In the second run (semi‐nested PCR), 1 μL from the first run PCR products was added into reaction mix containing Taq DNA polymerase (Qiagen). For group A and B viruses, AG655 and BG517 forward primers, respectively, and F164 reverse primers were used as described by Sullender et al.28 The length of amplified DNA products for A, B, and BA viruses were 450 bp, 585 bp, and 645 bp, respectively.

Before nucleotide sequencing, the amplified PCR products were purified using the Qiagen PCR purification kit or Qiagen Qiaquick Gel Extraction Kit (Qiagen) from agarose gel. Cycle sequencing was performed using a BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) in both directions. Sequence data were generated on the ABI 3130 Genetic Analyzer (Applied Biosystems). Sequences were manually edited, and multiple sequence alignments were generated using the Clustal X 2.012 version.29 Identical sequences were identified with DAMBE software, version 5.3.10.30

Alignment of sequences from our study with reference sequences was performed on the ClustalW programme with MEGA6 software.31 For HRSVA and B, the phylogenetic tree of HVR2 was constructed by the neighbour joining method using the substitution model nucleotide P‐distance, complete deletion for gap or missing data treatment and 1000 replicates of bootstrap probabilities within MEGA 6.

The partial nucleotide sequences of G gene, belonging to the six epidemic seasons, were deposited at GenBank under accession numbers KM266805‐KM267027 and MF177907‐MF177952. The HRSV reference sequences used in this study are shown in Supplement 1A,B.

Statistical analysis

SPSS v12.0 software (SPSS, Inc., Chicago, IL) was used for statistical analysis. The significance of differences in the rates among various groups was tested with the Student's t‐test and Fisher's exact test. Statistical significance was considered for a two‐sided value of P < 0.05.

RESULTS

Age distribution and seasonality

From 2009 to 2015, a total of 5166 samples were processed with real‐time PCR, and 1079 were HRSV positive. HRSV positivity according to years are given in Figure 1. Approximately 20% of the total samples was HRSV positive. The major populations afflicted by HRSV infections were infants (35.89%) and children between 13 and 24 months of age (46.95%). The occurrence of HRSV in the various age groups was as follows: 7.68% for 3‐5 years of age, 3.2% for 6‐13 years of age, 2.15% for 14‐27 years of age, 3.37% for 28‐60 years of age, and 0.65% for older than 60 years of age (Figure 2). Although HRSV activity starts at approximately the 47th week, the virus activity started at the 42nd week in the 2010‐2011 season. Usually, the epidemiological peak is reached at around week 52 in the winter (Figure 3). Both HRSVA and HRSVB co‐circulated in six epidemic seasons. HRSVA was dominant in the all epidemic seasons except 2011‐2012 season in which HRSVB was dominant genotype (Table 1).

Figure 1

HRSV positivity year‐based (2009‐2015)

Figure 2

Age distribution of HRSV‐positive samples

Figure 3

Year‐based (2009‐2015) weekly distribution of HRSV infections in Turkey

Table 1

HRSV‐positive samples and subgroup distribution in Turkey from 2009 to 2015

	No. RSV‐positive samples	No. samples typed	HRSVA	HRSVB
2009/2010	81	42	35	7
2010/2011	118	98	87	11
2011/2012	116	85	15	70
2012/2013	282	139	122	17
2013/2014	231	150	130	20
2014/2015	251	100	90	10

Clinical characteristics

Clinical data were available for 69 HRSV‐positive patients who were monitored due to ALRTI, especially bronchiolitis. Of 69 patients, 40 (58%) were male, and the median age was 4.9 months (10 days‐180 months) (IQR 25‐75%). Also considering the distribution of HRSV positivity according to gender, it was seen that males were more commonly infected than females (Figure 4).

Figure 4

HRSV positivity according to the gender

Thirteen patients had mixed viral infections: six patients had rhinovirus, one patient had adenovirus, two patients had human bocavirus, two patients had influenza B virus, one patient had parainfluenza virus 3 and one patient had parainfluenza virus 2+rhinovirus. Forty‐five (65.2%) of the patients were positive for HRSVA, and 24 (34.8%) were positive for HRSVB. Of the 45 HRSVA‐positive patients, 31 had NA1, and 14 had ON1, while patients with HRSVB had BA9. The most common symptoms observed in both subgroups included coughing (97.8% HRSVA and 100% HRSVB, respectively), wheezing (77.8% and 75%), rhinorrhoea (37.8% and 45.8%), and fever (37.8% and 29.2%). Patients in both groups had similar length of hospitalization, laboratory findings, and treatment (Table 2).

Table 2

Clinical symptoms and characteristics observed for HRSVA and B from the patients at the Pediatric Allergy and Immunology Clinic, Department of Children's Health Pediatric Hematology Oncology Training and Research Hospital, Ankara

	HRSVA (n = 45)	HRSVB (n = 24)	P
Age, month	5.6 (2.3‐10.6)	2.8 (1.8‐9.6)	0.20
Male gender, n (%)	27 (60.0)	13 (54.2)	0.64
Premature birth, n (%)	6 (13.3)	6 (25.0)	0.19
Clinical symptoms, n (%)
Fever	17 (37.8)	7 (29.2)	0.42
Coughing	44 (97.8)	24 (100)	0.65
Wheezing	35 (77.8)	18 (75.0)	0.79
Rhinorrhea	17 (37.8)	11 (45.8)	0.61
Treatment, n (%)
Oxygen	37 (82.2)	21 (87.5)	0.42
Antibiotics	36 (80.0)	21 (87.5)	0.33
Systemic steroid	23 (51.1)	15 (62.5)	0.36
Salbutamol	44 (97.8)	21 (87.5)	0.11
Duration of response to treatment, day	5 (3.75‐6.0)	5 (5.0‐8.0)	0.1
Length of hospitalization, day*	9 (7.0‐9.0)	8 (7.0‐11.0)	0.65
Laboratory findings, n (%)
Leukocytosis	5 (11.1)	1 (4.2)	0.33
Eosinophilia	1 (2.2)	1 (4.2)	0.59
CRP	17 (37.8)	12 (50.0)	0.25
Pathological findings in chest X‐ray	29 (64.4)	20 (83.3)	0.062

Phylogenetic analysis

Partial sequence analysis of the HVR2 region of G gene of 479 HRSVA and 135 HRSVB was aligned with representative sequences of identified reference genotypes. Phylogenetic trees are shown in Figures 5A and 5B. According to novel genotype criteria, proposed by Venter et al8 with grouping in the same phylogenetic cluster, a bootstrap value of 70‐100% and P distance values smaller than 0.07, the TR‐A genotype, belonging to the HRSVA genotype, and TR‐BA1 and TR‐BA2 genotypes, belonging to HRSVB genotype, were identified in our country between 2009 and 2013.8 The bootstrap values were 83% for TR‐A, 76% for TR‐BA1, and 75% for TR‐BA2 in clusters. The P‐distance values of the clusters were 0.01.

Figure 5

Unrooted phylogenetic trees of HRSVA (A) and HRSVB (B) strains from the Turkey. The phylogenetic tree of HVR2 was constructed by the neighbour‐joining method using the substitution model nucleotide p‐distance, complete deletion for gap or missing data treatment, and 1000 replicates of bootstrap probabilities within Mega 6

Cui et al13 described a novel genotype from China, named CB1. The CB1 genotype was defined by K224E, R214I, and V251M amino acid substitutions.13 Auksornkitti et al20 described a novel genotype from Thailand, called THB, with amino acid substitutions R214I, D215V, I227A, K234N, V251M, K258N, and E272D.20 Ren et al22 reported the GB5 genotype from China with the same amino acid substitution. The CB1, GB5, and THB strains clustered in the same branch with a bootstrap value of 100% (Figure 5B). Neither, Ren et al22 nor Auksornkitti et al20 include CB1 genotype in their studies. It is obvious that three genotypes are identical genotypes.

Amino acid analysis

Partial sequences of the HVR2 region of the G gene of 479 HRSVA and 135 HRSVB were successfully obtained. The mean nucleotide diversity within HRSVA is 0.062 and HRSVB is 0.058.

Four putative N‐glycosylation sites, N237, N251, N273, and N294, were identified in HVR2 of HRSVA (in relation to the HRSVA2 prototype strain). N237 was detected in only one NA strain and in the majority of the ON1 strains, and N251 and N273 were detected in most of the NA1 strains. N294 remained conserved in all NA1 and ON1 strains. In strain TR‐A, there were N251 and N294 N‐glycosylation sites. In the HVR2 region of Turkey strains, there were 32 Thr and Ser residues in the NA1 genotype, 34 Thr and Ser in TR‐A, and 42 Thr and Ser in ON1 that could be O‐glycosylation sites.

N‐glycosylation sites at the N276 and N290 positions were identified as conserved spots in the HRSVB strains from Turkey. In BA9 strains with nine successive nucleotide insertions, only N290 N‐glycosylation site was conserved. In TR‐BA1, there were N253, N273, N276, and N290 (relating to HRSVBA prototype strain) N‐glycosylation sites. In TR‐BA2, N276 and N290 N‐glycosylation regions were identified. In the HRV2 region of the HRSVB Turkey strains, there were approximately 48 Thr and Ser residues.

The predicted lengths of NA1 and TR‐A are 297 and 321 amino acids, respectively, in strain ON1. In a novel subgenotype, TR‐A, detected genotype‐specific substitutions include E232G, I243T, Q261P, S277P, and Y285H. Furthermore, the S250T substitution was observed in strains GA3, GA6, and GA7. E232G was reported as unique for ON1. Two other substitutions, T253K and S314L, were conserved in Turkey strains (Supplement 2).

The lengths of the G protein were 312, 315, and 319 amino acids in Turkey BA9 strains. In four Turkey BA9 strains, after the 693 amino acid position, insertion of nine successive nucleotides (ACAAAAAAA) coding Thr (T), Lys (K), Lys (K) (in relation to HRSV‐B prototype strain) was detected. Amino acid substitutions reported in BA9 were conserved in most of the Turkey strains (P219L, M222T, L237P, P258S, etc.). Novel genotype TR‐BA1 genotype‐specific substitutions were V251I, D253N, D273N, H287Y, and E292G and novel genotype TR‐BA2 genotype‐specific substitutions were E226D, I229T, T270I, and E292G (relating to the HRSVBA prototype strain) (Supplement 3).

The distribution of genotypes, belonging to each HRSV season, is given in Table 3. NA1 was observed as the dominant genotype during two seasons, BA9 was dominant in the 2011‐2012 season and ON1 was dominant in the 2012‐2013 season. The TR‐A sub‐genotype was detected during four successive epidemic seasons; TR‐BA1 and TR‐BA2 were detected in the 2010‐2011 and 2011‐2012 epidemic seasons.

Table 3

HRSV genotype distributions in Turkey from 2009 to 2015

	Genotypes
	NA1	ON1	BA9	BA10	TR‐A	TR‐BA
2009/2010	34	0	7	0	1	0
2010/2011	80	0	9	1	7	1
2011/2012	7	7	56	6	1	8
2012/2013	8	103	17	0	11	0
2013/2014	10	120	20	0	0	0
2014/2015	0	90	10	0	0	0

DISCUSSION

HRSV is the most common viral pathogen causing lower respiratory tract infections among infants and young children.32 The incidence of HRSV‐associated bronchiolitis and pneumonia was highest in 2‐ to 3‐month‐old infants; with increasing age, the incidence decreased.33 In this study, bronchiolitis was common for both subtypes according to clinical analysis of 69 pediatric patients. A higher susceptibility of HRSV infection occurred at the age of 0‐2 years. Most infected population were infants (35.4%) and children between 13 and 24 months (47.7%).

The HRSV prevalence decreased to 7.4% in children between 3 and 5 years old. The prevalence of HRSV decreased in children older than 1 year of age.34, 35 Lower immune response in the early years will develop soonest when children are re‐infected for the second or third time.36 There were no significant differences in the clinical symptoms, similar length of hospitalization, laboratory findings, and treatment observed between the HRSVA and HRSVB groups or in the coinfections in this study.

Although any clinical differences between the HRSV genotypes were not reported.11, 13, 18, 37, 38, 39 Baek et al11 noted that HRSVB infections had a relatively greater impact on the disease severity. Coinfection had no obvious impact on the disease severity.11 In our study, males were more commonly infected than females, as mentioned before.11, 39, 40, 41

RSV causes epidemics in the beginning of winter and spring in temperate climates, but its activity may differ locally.42, 43, 44, 45 HRSV infections appear at the end of November and peak in January in Turkey. Usually, an epidemiological peak is reached at around week 52 in the winter.

Considering the genotype analyses, both genotypes co‐circulated in six studied epidemic seasons. HRSV group dominance has been reported worldwide; both genotypes should be co‐circulated and dominant genotype should be variable. In many studies, HRSVA viruses are more frequently detected.13, 18, 39 In Turkey, HRSVA was dominant in two successive epidemic seasons. HRSVA was replaced by HRSVB in the 2011‐2012 season in which BA9 was dominant. Then, in the next season, HRSVB was replaced by HRSVA in the 2012‐2013 epidemic season. This replacement between subgroups affected by various factors, such as adaptation to host, antigenic variants, and social immunity.11, 41, 46

Although TR‐A, TR‐BA1, and TR‐BA2 were considered emerging subgenotypes in this study, they were detected in very few samples. Until 2013, the dominant genotype of HRSVA was NA1 in Turkey. Turkey NA1 strains lost the N237D N‐glycosylation site, except one strain. In TR‐A, there were genotype‐specific amino acid substitutions (I243T, S277P, S250T, and Y284H). Two N‐glycosylation sites were detected. The I243T substitution was a potential O‐glycosylation site. ON1 was identified in the 2011‐2012 epidemic season in Turkey and became the dominant subgenotype in 2013. It is required to monitor whether these amino acid substitutions would affect the spread of virus.

Only BA9 of HRSVB was identified in Turkey. No other HRSV B genotype has been detected in our country since 2009. BA9 was dominant in Japan, Korea. and China during the same period.11, 13, 17 In TR‐BA1, there were four N‐glycosylation sites. TR‐BA1 and TR‐BA2 were detected in 2011‐2012 season, which is when BA9 was dominant. Such prevalence of BA9 in society would provide social immunity and cause changes in the viral structure. Similarly, TR‐A was detected in years when HRSVA was dominant.

The length of the G protein of HRSVB is variable due to changes in the amino acid substitution, insertion and stop codon. The polymorphism of the G protein of HRSVB is more than that of HRSVA. In four BA9 strains, ANK/118/2012, ANK/131/2012, ANK/138/2013, and ANK/507/2013, nine nucleotide insertions after the 693rd nucleotide were present. Insertion in the same place was determined by Cui et al13 in CB1 strains. The association of changes in the G protein length with antigenic variation has been reported in many studies. The antibody response induced by HRSVA may provide cross protection against HRSVB such that HRSVB undergoes more extensive antigenic variation to evade herd immunity.2, 47

Some recently reported novel variants, such as BAC and CB‐A within HRSV subgroups, seem to be circulating domestically.11, 13 TR‐A, TR‐BA1, and TR‐BA2 should also remain domestically. As proposed by Reiche and Schweiger,35 these results indicate that HRSV strains may develop independently. It should be monitored whether these novel genotypes were temporary or will establish a stable linage in the population. Additionally, it is important to follow‐up on the mutations that the virus developed against herd immunity.

Diversity in the HVR2 region of the G protein plays a significant role in HRSV pathogenesis by facilitating evasion of the immune response.8, 48 Our results indicate that the ectodomain G protein of HRSV has both N‐linked and O‐linked glycosylation with O‐linked sugars contributing to the major part of glycosylation.2, 7, 11

The THB and GB5 genotypes should be classified into the CB1 genotype. We strongly recommend using the CB1 genotype due to the priority rule. A similar issue was detected with the CB‐B and BA12 genotypes.11, 19 CB‐B and BA12 strains clustered in the same branch with a bootstrap value of 65% (Figure 5B). Additionally, the BA11 genotype was clustered in the BA9 genotype in our study.23

The present study provides additional information on the molecular epidemiology of HRSV in Turkey and presents data for comparative analysis with HRSV stains circulating in local communities and other regions worldwide. These findings may improve our understanding of HRSV evolution and the potential development of an HRSV vaccine.

Additional Supporting Information may be found online in the supporting information tab for this article.

Supporting Tables S1.

Click here for additional data file.

Supporting Data S1.

Click here for additional data file.