Journal of Medical Virology 86:1971–1982 (2014)

Molecular Epidemiology of Human Respiratory Syncytial Virus Over Three Consecutive Seasons in Latvia Reinis Balmaks,1,2 Irina Ribakova,1,2 Dace Gardovska,1,2 and Andris Kazaks3* 1

Department of Pediatrics, Riga Stradins University, Riga, Latvia Children’s Clinical University Hospital, Riga, Latvia 3 Latvian Biomedical Research and Study Centre, Riga, Latvia 2

Lower respiratory tract infections caused by the human respiratory syncytial virus (HRSV) represent an immense burden of the disease, especially in young children. This study aimed to investigate the evolutionary history of HRSV strains isolated in the Children’s Clinical University Hospital (Riga, Latvia) over three consecutive HRSV seasons. Of 207 samples from children hospitalized with lower respiratory tract infections, 88 (42.5%) tested positive for HRSV by RT-PCR. The seasonal activity started and peaked later than the average for the Northern hemisphere. Patients with HRSV lower respiratory tract infection were significantly younger than patients not infected with HRSV. HRSV-A viruses predominated for two consecutive seasons and were followed by an HRSV-B dominant season. Phylogenetic analysis based on glycoprotein G gene partial sequences revealed that viruses of both groups belonged to the worldwide dominant genotypes NA1 (HRSVA) and BA-IV (HRSV-B). High diversity of this gene was driven only partially by selection pressure, as only two positively selected sites were identified in each group. Two of the HRSV-A isolates in this study contained a 72nt duplication in the C-terminal end of the G gene (genotype ON1) that was first described in Canada in the 2010–2011 season. Initial spatial and temporal dynamics of this novel genotype were reconstructed by discrete phylogeographic analysis. Fifteen years after acquiring comparable 60-nt duplication in the G gene, genotype BA lineages have replaced all other HRSV-B strains. However, the population size of genotype ON1 plateaued soon and even decreased slightly before the beginning of the 2012–2013 season. J. Med. Virol. 86:1971–1982, 2014. # 2013 Wiley Periodicals, Inc. C 2013 WILEY PERIODICALS, INC. 

KEY WORDS:

HRSV genotypes; glycoprotein G; phylogenetic analysis; phylogeographic analysis

INTRODUCTION The global burden of lower respiratory tract infection caused by human respiratory syncytial virus (HRSV) is immense. During annual epidemics it is a major reason for pediatric hospitalizations worldwide and in the developing countries it is also one of the leading causes of childhood mortality [Hall et al., 2009; Lozano et al., 2012]. A substantive disease burden is also associated with HRSV in vulnerable adult populations [Falsey et al., 2005]. Although there is an obvious medical need, no licensed vaccine is available yet [Anderson et al., 2013]. HRSV is a member of the Pneumovirus genus that is classified within the subfamily Pneumovirinae of the family Paramyxoviridae. Accordingly, it is a cytoplasmic, enveloped virus with linear, negative sense, ssRNA genome [Wang et al., 2012]. The viral RNA of HRSV is approximately 15.2 kB in size and encodes 11 viral proteins [Collins and Crowe, 2007]. Two surface glycoproteins, G and F, are antigenically significant because they induce neutralizing antibody Grant sponsor: The Study and Science Administration Grant No 09.1604 sponsored by Ministry of Education and Science of the Republic of Latvia; Grant sponsor: European Social Fund Grant No 2009/0147/1DP/1.1.2.1.2/09/IPIA/VIAA/009 The authors declare that they have no conflict of interest.
Correspondence to: Andris Kazaks, Latvian Biomedical Research and Study Centre, Ratsupites 1, Riga, LV-1067, Latvia. E-mail: [email protected] Accepted 29 October 2013 DOI 10.1002/jmv.23855 Published online 2 December 2013 in Wiley Online Library (wileyonlinelibrary.com).

1972

responses [Connors et al., 1991]. Based on the reaction with monoclonal antibodies, HRSV strains are separated into two major groups, HRSV-A and HRSV-B [Mufson et al., 1985], which are genetically divergent viruses that have evolved separately for about 350 years [Zlateva et al., 2005]. Viruses from both groups co-circulate simultaneously during epidemic seasons with alternating predominance. Typically, there is a cyclic pattern whereby several predominant HRSV-A seasons are followed by a single HRSV-B dominant season [Venter et al., 2001; Zlateva et al., 2007]. HRSV viruses also vary considerably within the groups, with several distinct genotypes in each group accounting for clusters of circulating strains. Several genotypes co-circulate in the same community and are replaced by new ones in successive seasons [Peret et al., 2000; Venter et al., 2001; Zlateva et al., 2007]. The new genotypes spread worldwide very quickly, and are related more temporally than geographically [Peret et al., 1998; Trento et al., 2010]. The most extensive differences are found in the gene encoding G protein, and the genotype classification based on partial sequencing of this gene is used widely in molecular epidemiologic studies of HRSV [Peret et al., 1998, 2000]. It has been confirmed by genome wide analysis that genotyping based G gene variability represents overall virus variability [Tan et al., 2012]. In the reference strain A2, the G protein is a type II transmembrane glycoprotein 298 amino acids (aa) in length. Moving in the direction from N- to Cterminal end, there are three structural domains: cytoplasmic (aa 1–37), transmembrane (aa 38–63), and ectodomain (from aa 64). The large ectodomain consists of two mucin-like hypervariable regions (HVR1 & 2) separated by a segment of 13 highly conserved aa (positions 164–176 in strain A2) [Johnson et al., 1987]. Sequence diversity in the two hypervariable regions is among the most extensive found in human viruses. Significant alterations such as aa substitutions, deletions, insertions, frame-shift mutations, large duplications, and premature termination codons have been reported, representing both immune-driven selection and structural plasticity of the protein [Trento et al., 2003; Zlateva et al., 2005; Collins and Melero, 2011]. Thus, the resulting protein length can vary from 254 to 321 aa [Zlateva et al., 2005; Eshaghi et al., 2012]. The mucin-like regions, like eukaryotic mucins, are rich with serine, threonine and proline residues, and are glycosylated heavily with O- and N-linked carbohydrate side chains. The number and positions of glycosylation sites are conserved poorly among the strains, also contributing to their antigenic differences [Johnson et al., 1987; Palomo et al., 1991; Martinez et al., 1997]. Local molecular surveillance is important to identify and predict new virulence factors, select vaccine candidate strains and develop public health intervention strategies based on phylogeographic analysis. The aim of this study was to investigate the history J. Med. Virol. DOI 10.1002/jmv

Balmaks et al.

of and the processes driving the evolution of HRSV strains isolated in the Children’s Clinical University Hospital (Riga, Latvia) over three consecutive seasons. METHODS Patients, Isolates, HRSV Detection and Differentiation The study protocol was approved by the Ethics Committee of Riga Stradins University, and written informed consent was obtained from the parents of all enrolled children. Inclusion criteria, sample collection, HRSV detection, and group differentiation have been described previously in detail [Balmaks et al., 2011, 2013]. In brief, nasopharyngeal aspirates were collected from previously healthy 2–24-monthold children hospitalized in Children’s Clinical University Hospital in Riga, Latvia with lower respiratory tract infection over three consecutive seasons (2009–2012). Clinical samples were frozen and stored at 70˚C. Total RNA was extracted directly from nasopharyngeal aspirates and used as a template for cDNA synthesis by reverse transcription using primer F164 [Sullender et al., 1993]. HRSV was detected by PCR-amplification of a conserved fragment in the non-coding sequence between the P and M genes. HRSV-positive samples were differentiated into groups A and B by group-specific PCR targeting the HVR2 segment of the G gene. The amplified fragments from the isolates of seasons 2010–2011 and 2011–2012 were purified and sequenced by groupspecific forward primers 50 -CATATGCAGCAACAATCCAAC-30 and 50 -CCAATCCACACAAATTCAGC-30 for groups A and B, respectively and by a cross-reactive reverse primer 50 -CTCCATTGTTATTTGCCCCAG-30 . The reaction conditions were described previously [Balmaks et al., 2011]. The sequences were processed using Vector NTI Advance 10 (Invitrogen, Carlsbad, CA) software package. A 336-nt-long HRSV-A sequence (corresponding to codon positions 187–299 of reference strain A2; GenBank accession number M11486) and a 516-nt-long HRSV-B sequence (corresponding to codon positions 140–293 of reference strain B1; GenBank acc. no. AF013254) were deposited in GenBank under acc. nos. KF030137–KF030185. Full length G gene sequences from season 2009–2010 have been published before (GenBank acc. nos. JF979145–JF979157) [Balmaks et al., 2013]. A nomenclature was adapted for Latvian isolates where LV indicates Latvia followed by the laboratory log number and the isolation year. Phylogenetic and Adaptive Evolutionary Analysis The alignments of nucleotide (nt) and deduced aa sequences were prepared via the EMBL-EBI ClustalW2 software [Larkin et al., 2007]. The phylogenetic trees were constructed using the neighbor-joining algorithm and genetic distances (the number of nt

Three-Year Epidemiology of HRSV in Latvia

and aa substitutions per site from averaging over all sequence pairs) were calculated under the best-fit substitution models in MEGA5.1 software [Tamura et al., 2011]. Bootstrapping with 1,000 replicates was performed for each analysis to evaluate confidence estimates. Nucleotide diversity (the average number of nt differences per site; p) was calculated using DNAsp v5 software [Librado and Rozas, 2009]. The Oglycosylation sites were determined using NetOGlyc 3.1 server neural network predictions; a G-score >0.5 was used as the cut-off [Julenius et al., 2005]. Acceptance of the N-linked oligosaccharides of the consensus motif NXS/T, where X is any aa except proline [Hart, 1992], was predicted when glycosylation potential was >0.5 in NetNGlyc 1.0 server [Gupta et al., 2004]. Positively selected sites were identified by estimating site-specific non-synonymous (dN) to synonymous (dS) substitution rate ratios (dN/dS, or v) with five different algorithms available on the Datamonkey web server [Delport et al., 2010]: single likelihood ancestor counting (SLAC), fixed effects likelihood (FEL), internal fixed effects likelihood (IFEL), random effects likelihood (REL), and mixed effects model of evolution (MEME). The best-fit substitution model was determined using MEGA5.1 software. The site was considered under positive selection (dN/dS > 1) when two or more methods reached agreement with statistical significance (P < 0.1 or Bayes factor >20). The mean dN/dS ratio was estimated using the SLAC algorithm. Analogous gene fragments, that is, codon positions 187 to Stop of reference strains A2 and B1, were used, whenever parameters of the two groups were compared. Evolutionary Rate, Population Dynamics and Phylogeographic Analysis Nucleotide substitution rate per site, the time of most recent common ancestor (tMRCA), changes in the population size and discrete phylogeographic

1973

analysis were estimated from time-stamped sequences with location trait using Bayesian Markov chain Monte Carlo (MCMC) method in BEAST v2.0.2 software package [Drummond et al., 2012]. The fittest nucleotide substitution model was determined using MEGA5.1. The dataset was then analyzed assuming relaxed (uncorrelated lognormal) molecular clock and using Bayesian skyline demographic model as coalescent prior. MCMC chain length was chosen to reach effective sample size >200. The results obtained from MCMC analysis were assessed using Tracer v1.5 and the maximum clade credibility (MCC) tree was inferred using TreeAnnotator v2.0.2. The MCC tree was visualized using FigTree v1.4.0 software. RESULTS HRSV Detection and Epidemics From July 1, 2009 to June 30, 2012 a total of 207 specimens were collected from hospitalized patients. Most of the patients were admitted to the general pediatric floor and 12 (6%) were admitted to the pediatric intensive care unit. The proportion of HRSV-positive samples in each season ranged from 33% to 57% (Fig. 1). Every year HRSV outbreaks started in winter and extended in spring. Overall, HRSV epidemic activity peaked at week 10 and ranged from weeks 46 to 21. Based on medical chart review, 5 (6%) of patients with HRSV infection had another viral co-infection detected by off-study investigations—2 influenza A, 2 influenza B, and 1 human metapneumovirus. The majority of HRSV-positive individuals were 0–6 months old (58%) and the number of HRSV-positive individuals decreased with increasing age. Patients with HRSV infection were significantly younger than those who tested negative (median age 6 months vs. 10 months; P < 0.001). Six patients were enrolled twice in the study, but none of them had HRSV re-infection.

Fig. 1. Number of total samples tested (curve) and number of HRSV-A (white bars) and B (black bars) positive samples per month among hospitalized patients in Latvia over three consecutive seasons. Months are indicated by the initial letter starting from July, 2009.

J. Med. Virol. DOI 10.1002/jmv

1974

Balmaks et al. TABLE I. Group Distributions Over HRSV Seasons No. (%) within HRSV-positive

Epidemic season

Samples tested

2009–2010 2010–2011 2011–2012 Total

No. (%) of HRSV-positive

69 67 71 207

23 38 27 88

Group Prevalence and Phylogenetic Analysis All isolates were typed according to HRSV groups and there were no cases of co-infection with the two groups. HRSV-A and B viruses co-circulated in all three seasons, however, HRSV-A predominated in seasons 2009–2010 and 2010–2011, while HRSV-B predominated in season 2011–2012 (Table I). The HVR2 segment of the G gene was sequenced for all HRSV-positive samples. By comparing nt sequences, 53 HRSV-A and 35 HRSV-B isolates grouped into 29 and 23 different strains, respectively. Representative gene sequences of strains retrieved in this study were aligned and included in phylogenetic analysis together with sequences of reference strains obtained from GenBank representing 11 HRSV-A genotypes (GA1– 7, SAA1, NA1 and 2, and ON1) and 13 HRSV-B genotypes (GB1–4, SAB1–3, BA-I–VI). Two HRSV-A genotypes were identified, with 28 strains (51 isolates) clustering in genotype NA1 (Fig. 2a). One strain representing two identical isolates from March, 2012 clustered in the novel, NA1-related, genotype ON1, characterized by a 72-nt duplication in the HVR2 [Eshaghi et al., 2012]. All HRSV-B strains had a 60-nt duplication in the G gene, characteristic of the BA genotypes, and clustered within the clade BAIV (Fig. 2b). Eight HRSV-A and one HRSV-B strains were found to be identical to previously elsewhere isolated sequences through NCBI-BLAST nucleotide query (Table II) [Zhang et al., 2000]. Molecular Analysis of Genotype NA1 The HVR2 sequences of the G gene of Latvian strains of genotype NA1 were related closely, with 2.1% (0.4) and 4.1% (0.9) divergence at nt and aa levels, respectively. The nucleotide polymorphism at codon sites and resulting aa alignment are shown in Figure 3a. Overall, in the 112-aa long sequence analyzed, 63 aa matching those of strain A2 were conserved across all isolates. All HRSV-A strains had a premature termination codon (TGA) at aa position 298 and preserved the original termination codon (TAG) at position 299. These strains also had the following G protein mutations: P226 ! L, E233 ! K, L258 ! H/Y, M262 ! E/R, F265 ! L, S269 ! T, S280 ! Y, P289 ! S, S290 ! P/L, P292 ! S, P293 ! S, P296 ! T, and R297 ! K. Other common mutations (detected in >50% of strains) included L208 ! I, S222 ! P, N237 ! D, I244 ! R, J. Med. Virol. DOI 10.1002/jmv

(33) (57) (38) (42.5)

HRSV-A 15 28 10 53

(65) (74) (37) (60)

HRSV-B 8 10 17 35

(35) (26) (63) (40)

N273 ! Y/K, P274 ! L, and P286 ! L. The mean divergence between Latvian strains and the NA1 reference strain NG-016-04 (GenBank acc. no. AB470478) was 2.3% (0.4) and 4.2% (0.9) at nt and aa levels, respectively, which was similar to the divergence within Latvian strains. Using NetOGlyc software, 40 serine and threonine residues of the consensus sequence of Latvian strains were predicted to be potentially O-glycosylated (range 36–43 among individual strains). In the same sequence strain A2 is predicted to have 39 O-glycosylation sites, with 36 of them being in common with the Latvian consensus sequence. There are two potential N-glycosylation motifs (NXS/T) in the HVR2 segment of strain A2 (N237 and N251), and two additional motifs were identified among Latvian strains (N273 and N294). NetNGlyc software predicted only N251 and N294 glycosylation of the Latvian consensus sequence. Among individual strains, N-glycosylation sites ranged from 0 to 3. The NLS tripeptide at position 273 is followed by proline which makes glycosylation of this site unlikely [Gupta et al., 2004], and glycosylation of NTT at position 294 appeared to depend on the P286 ! L mutation. None of the sites was conserved among all isolates. To identify positively selected sites, five different methods were used. In total, 28 sites showed elevated rates of non-synonymous substitutions (dN) and the mean non-synonymous/synonymous substitution rate ratio (dN/dS) was 1.12, suggesting presence of selection pressure. However, only two positively selected sites, M262 and N273, were identified based on the consensus of two or more of the methods (Table III, Fig. 3). The N273 substitution to tyrosine led to loss of the NXS/T motif. Molecular Analysis of Genotype BA-IV Latvian strains of the BA-IV genotype were more diverse than NA1 strains, with calculated mean distances of 3.3% (0.4) and 4.2% (0.7) at nt and aa levels, respectively. This is also evident by nucleotide polymorphism analysis (Fig. 3b); 23 codons had p values >0.1. In addition to the characteristic 20-aa in-frame duplication in the HVR2 segment, all Latvian BA-IV strains had deletion P159–K160. This 2-aa deletion offset all the codon positions which in this article are given according to the reference strains. Following this deletion there was an absolutely conserved 28-aa region (aa positions 161–188)

Three-Year Epidemiology of HRSV in Latvia

1975

Fig. 2. Neighbor-joining phylogenetic trees for HRSV-A (A) and B (B). The dendrograms were calculated based on alignments of non-redundant sequences of previously published (^) and Latvian strains isolated in this study (^) together with reference strains with assigned genotypes. TN93 substitution model and discrete gamma distribution with four rate categories were used for both groups. The genotypes are labeled as follows: GA1–7 and GB1–4 according to [Peret et al., 2000],

SAA1 and SAB1–3 according to [Venter et al., 2001], NA1 and 2 according to [Shobugawa et al., 2009], ON1 according to [Eshaghi et al., 2012], and BA-I–VI according to [Trento et al., 2006]. Taxa description include: strain name in bold, GenBank acc. no. in brackets and number of identical isolates in italics. Bootstrap values of 70 are shown at nodes; both trees are drawn in the same scale and the bar denotes nucleotide substitutions per site.

recognized in all Latvian strains. Two alternative terminal codon positions were observed either at position 293 (TAA) or at the original position 300 (TAG) of the reference strain B1 (299-aa-long). Thus, the Latvian HRSV-B viruses produced either a 310or 317-aa-long G protein. The 317-aa-long strains emerged only in season 2011–2012. Alignment of protein sequences together with the originally described BA strain BA/3833/99 is shown in Figure 3b. Overall, in the 171-aa-long sequence analyzed, 123 aa of strain BA/3833/99 were conserved across all isolates. All Latvian BA-IV strains had K218 ! T, L223 ! P, and S247 ! P mutations. In addition, the following mutations were very common (>50%): L219 ! P, T270 ! I/F, V271 ! A, H287 ! Y, and Q313 ! Stop. The mean divergences between Latvian BA strains and the reference strain BA/3833/99 were 4.2% (0.6) and 5.2% (0.9) at nt and aa levels,

respectively. NetOGlyc software predicted 44 O-glycosylation sites in the consensus sequence of Latvian strains (range 36–50) and 41 in strain BA/3833/99. Sixteen of the predicted O-glycosylation sites were in common between consensus sequences of HRSV-A and B. There are two potential N-glycosylation motifs in the HVR2 segment of strain BA/3833/99 (N296 and N310). N296 was conserved and predicted to be glycosylated in all strains, while N310 was predicted to be glycosylated only in 317-aa-long sequences. Strain LV/057/11 also had three additional NXS/T motifs (N212, N230, and N253), two of which (N212 and N230) were predicted to be glycosylated. Thirtytwo sites showed elevated rates of dN and the overall dN/dS ratio was 0.45. Two positively selected sites, L219 and T270, were found by two-method agreement. The latter codon is situated in the duplicated region (Table III, Fig. 3). J. Med. Virol. DOI 10.1002/jmv

1976

Balmaks et al. TABLE II. Previously Published Strains Identical to Isolates From Latvia in 2009–2012 Latvia

World Season

Straina

Season

Country

GenBank

LV/009/11 LV/011/11

2010–2011 2011–2012 2010–2011 2010–2011

LV/027/12 LV/036/12

2011–2012 2011–2012

LV/043/12

2011–2012

LV/046/12 Genotype ON1 LV/029/12

2011–2012

NL20752227/07-08 Riyadh 42/2008 08-047045 A/RJ/64/2006 A/RJ/288/2007 A/RS/5357/2009 A/WI/629-22/07 R9061/07-08 R14944/08-09 18227AN/2012 18294RM/2012 HD/70178486/11 ON69-0310A ON187-0111A A/WI/629-Q0284/10 07408AN/2012 HD/70189233/12 16158AN/2012

2007–2008 2008–2009 2008–2009 2006–2007 2007–2008 2009–2010 2007–2008 2007–2008 2008–2009 2011–2012 2011–2012 2011–2012 2009–2010 2010–2011 2009–2010 2011–2012 2011–2012 2011–2012

The Netherlands Saudi Arabia The Netherlands Brazil Brazil Brazil USA UK UK Italy Italy Germany Canada Canada USA Italy Germany Italy

HQ731761 JX131640 JX015483 JX182800 JX182830 JX182837 JF920049 HQ731741 HQ731736 JX988459 JX988494 JX967562 JN257697 JN257701 JF920053 JX988460 JX967566 JX988467

ON138-0111A 0284RM/2012 WUE/9583/12 RSVA/GN435/11 RSV/Yokohama.JPN/P6587/2012

2010–2011 2011–2012 2011–2012 2011–2012 2011–2012

Canada Italy Germany South Korea Japan

JN257694 JX988452 JX912363 JX627336 AB761611

BE/8844/09

2009–2010

Belgium

JX645922

Strain Genotype NA1 LV/056/10

2011–2012

Genotype BA-IV LV/010/11 2010–2011 a

Strains with 100% maximal identity in NCBI-BLAST query. Only one strain per country per season is listed.

Phylodynamics and Phylogeography of the Novel Genotype ON1 Two identical HRSV-A isolates (LV/029/12) belonged to the novel genotype ON1 and were identical to the originally isolated sequence. Genotype ON1 is characterized by 72-nt duplication in the HVR2 that leads to indel mutation E284 ! G and duplication of Q261–S283, elongating G protein by 24 aa [Eshaghi et al., 2012]. Because of the premature termination codon like in NA1 strains the total protein length is 321 aa. The duplicated region was predicted to contain seven additional O-glycosylation sites, while it interfered with N-glycosylation of N318(294) (Fig. 3). Isolates of this study came from two unrelated patients in March 2012. The patients were admitted to the hospital at different times and were located on different floors, that is, were unlikely to have an intra-hospital transmission. Since the discovery in Canada in 2011 this novel genotype has been isolated in at least nine different countries in North America, Europe, Africa, Asia, and Oceania with the total of 42 sequences available in the GenBank [Eshaghi et al., 2012; Lee et al., 2012; Cui et al., 2013; Khor et al., 2013; Prifert et al., 2013; Valley-Omar et al., 2013]. Sixteen of these sequences were identical to the strain ON67-1210A and overall showed little diversity with the mean divergence of 0.8% (0.2) and 1.5% (0.5) at nt and aa levels, J. Med. Virol. DOI 10.1002/jmv

respectively. The evolutionary rate of the HVR2 was calculated to be 7.92  103 nt substitutions/site/year (highest probability density; HPD 95%: 2.97  103– 1.28  102) and the tMRCA was estimated August 2010 (HPD 95%: September 2009–January 2011). Discrete phylogeographic analysis depicted by the MCC tree demonstrates a hypothesis of global dissemination (Fig. 4). This analysis suggests spread from Canada to Italy and Germany, initially, and subsequent spread to the rest of the world (posterior value >0.7). Also, dissemination from Germany to Japan and from South Africa to China is supported well. However, the migration of this virus to Latvia through Malaysia has poor support (posterior value 0.03). The Bayesian skyline model indicates that the effective population size of genotype ON1 had slow expansion, reached plateau early in 2012 and decreased slightly before the beginning of the season 2012–2013 (Fig. 4). DISCUSSION This study of HRSV variability was conducted in the only tertiary level care pediatric hospital in Latvia; therefore, results obtained are likely to be representative of the entire country. 42.5% of children 2 years of age hospitalized with lower respiratory tract infection in Latvia were diagnosed with HRSV which is the same as the average

Three-Year Epidemiology of HRSV in Latvia

1977

Fig. 3. Polymorphisms of the coding sequences and the alignments of non-redundant amino acid sequences of the G protein of HRSV-A (A) and B (B) strains isolated in Latvia. The x-axis of the graph represents the codon and the corresponding aa positions and the y-axis represents sequence polymorphisms expressed as nucleotide diversity p. The alignments are shown relative to the sequences of reference strains A2 (GenBank acc. no. M11486) of HRSV-A and BA/3833/99 (GenBank acc. no. AY333362) of the HRSV-B

genotype BA. Consensus sequences were generated using Align X software in Vector NTI Advance package. Identical sequences are identified as dots. Asterisks indicate the positions of stop codons. The aa duplication regions are enclosed in open boxes. Predicted O-linked glycosylation sites are indicated by (5) for reference strains and (!) for Latvian consensus. NXS/T motifs that were predicted to be Nglycosylated are shaded gray. Positively selected sites are marked by (þ).

TABLE III. Sites Found to be Under Positive Selection

in Europe (42–45%) [Simoes and CarbonellEstrany, 2003]. A lower isolation rate in 2009–2010 [Balmaks et al., 2011] could be explained by the influenza A virus H1N1 pandemic that season [Nikiforova et al., 2011]. Peak seasonal activity in Latvia started consistently later in winter and lasted later in spring (February–April) compared to the average (December–February) of the Northern hemisphere [Bloom-Feshbach et al., 2013], but was similar to that observed in Russia [Tatochenko et al., 2010] and in biennial “late” seasons in Sweden [Reyes et al., 1997]. The reasons for geographic HRSV seasonality differences are not defined well, but are thought to be

Codon HRSV-A M262 N273 P274 HRSV-B L219 T270

SLAC

FEL

IFEL

REL

MEME

0.45 0.59 0.26

0.16 0.08 0.21

0.03 0.09 0.07

2.89 3.42 2.95

0.01 0.08 0.23

0.06 0.30

0.09 0.05

0.04 0.14

14.36 65.13

0.11 0.07

SLAC, FEL, IFEL, and MEME significance levels are given as Pvalues, while empirical Bayes factors are given for REL. Bold values are considered statistically significant (P < 0.1, BF > 20). Underlined sites meet the criteria for two or more methods.

J. Med. Virol. DOI 10.1002/jmv

1978

Fig. 4. Discrete phylogeographic analysis and demographic history of genotype ON1. Forty-three partial HRSV-A G gene sequences with 72-nt duplication retrieved from GenBank were aligned and used to infer maximum clade credibility (MCC) tree (Upper Panel) and the Bayesian skyline plot (BSP; Lower Panel). HKY substitution model and discrete gamma distribution with four rate categories were used. Fifty million MCMC samples were generated and 1,000 were discarded as burn-in to reach convergence. Branch colors and widths in the MCC tree represent location (see legend) and posterior support for the branch (values 0.7 are shown at nodes), respectively. The yaxis of the BSP represents the population size, which is equal to the product of the effective population size (Ne) and the generation length in years (t). The red solid line represents the median estimate and the area between the dotted lines represents the 95% HPD limits. Both graphs are drawn in the same time scale. Taxa names code for country of isolation, isolation number and time.

related to local climate and host behavior factors [Stensballe et al., 2003]. HRSV seasons were as long as 6 months, which correlates with the global median [Bloom-Feshbach et al., 2013]. A larger community study is necessary to strengthen these data, which are crucial for local timing of immunoprophylaxis with palivizumab for high-risk infants. Patients with HRSV infection were significantly younger and were mostly less than 6 months of age. Viruses of both HRSV groups co-circulated in each of the three seasons in Latvia that were studied. HRSV-A strains were isolated more frequently and dominated in the first two of the analyzed seasons. Different patterns of group circulation have been J. Med. Virol. DOI 10.1002/jmv

Balmaks et al.

described, but generally 1–3 successive predominant HRSV-A seasons are followed by a single season of HRSV-B dominance [Venter et al., 2001; Zlateva et al., 2007]. The reason for this alternating pattern is not entirely clear, but lower host susceptibility to homologous than heterologous strains has been observed and suggested to play a role [White et al., 2005]. Peak activity of both viruses occurred simultaneously and there was no correlation between epidemiologic week and the HRSV group (data not shown). For in-depth molecular analysis, a segment (HVR2) of the highly variable glycoprotein G gene was sequenced from all isolates from Latvia. Multiple identical sequences to strains isolated in various countries were detected, suggesting that outbreak strains have no geographic restrictions and can remain stable for several epidemiologic seasons. For example, strain LV/056/10, represented by nine identical isolates, was detected in Latvia for two seasons, but has been circulating globally since the 2007–2008 season (Table II, Fig. 2). Phylogenetic analysis revealed that strains circulating in Latvia during the study period belonged to three different genotypes. A single genotype in each group, NA1 in HRSV-A and BA-IV in HRSV-B, remained dominant in all three of the seasons. In the previous study in Latvia, NA1 genotype was not included the analysis [Balmaks et al., 2013]. This genotype was described first in Japan during the 2004–2005 season [Shobugawa et al., 2009], and was reported subsequently in numerous molecular epidemiologic studies as the dominant HRSV-A genotype [Arnott et al., 2011; Eshaghi et al., 2012; Forcic et al., 2012; Ohno et al., 2013]. The NA1 genotype is related closely to GA2 and is not distinguished universally [van Niekerk and Venter, 2011; Yoshida et al., 2012]. All HRSV-B sequences had a 60-nt duplication in the G gene, characteristic of BA (Buenos Aires) genotypes. Since 2005 BA-IV has been the most common HRSV-B genotype detected worldwide and has replaced all other lineages [Agrawal et al., 2009; Trento et al., 2010; van Niekerk and Venter, 2011]. Additional sublineages within BA-IV have been reported [Dapat et al., 2010], but inclusion of them in the analysis led to inconsistent results. In 2012, a year after its first description, another HRSV-A genotype, ON1, was detected. ON1 is characterized by a 72-nt in-frame duplication in the G gene. This strain emerged in Canada in the 2010– 2011 season, where it accounted for 10% of HRSV-A infections [Eshaghi et al., 2012]. A similar distribution was reported from Germany in the following season [Prifert et al., 2013]. In the 2011–2012 season in Latvia only ten HRSV-A viruses were detected, of which two (20%) were ON1. Unlike most proteins, the divergence of HRSV G protein was greater at the aa level than at the nt level, confirming previous observations [Johnson et al., 1987; Cane, 1997; Venter et al., 2001]. This

Three-Year Epidemiology of HRSV in Latvia

underlines the structural tolerance of aa substitutions in this gene. There was a significant G protein length polymorphism among the isolates, with predicted lengths varying from 297 to 321 aa, accounted for by the following: (i) all NA1 sequences had a premature termination codon and were 297 aa long; (ii) the ON1 strain had a 321-aa-long sequence, resulting from a 24-aa insertion; (iii) BA-IV lengths were 310 and 317 aa, resulting from the use of alternative termination codons. The use of different reversible termination codons, including Stop313 ! Q, in HRSV-B has been described before [Botosso et al., 2009]. Together with in-frame aa insertions and deletions (such as deletion of P159–K160 observed in all Latvian isolates) there are at least 11 different protein lengths in BA lineages [Zlateva et al., 2007; Trento et al., 2010]. In vitro, the reduction of G protein size by in-frame premature termination codons allows neutralization escape by monoclonal antibodies without affecting HRSV infectivity [Rueda et al., 1991]. Strain specific monoclonal antibodies react only with glycosylated G protein [Martinez et al., 1997], therefore the glycosylation properties were evaluated by software predictions. Because of substitutions directly in glycosylation sites or surrounding sequences, O- and N-glycosylation sites were conserved poorly among the strains of both groups, further diversifying strain phenotypes. Despite sequence diversity and length polymorphism among the strains, the overall proline/serine/threonine proportions remained stable (roughly 10, 10, and 25%, respectively), maintaining mucin-type characteristics [Julenius et al., 2005]. In-depth natural selection analysis was also performed. Although the overall dN/dS ratio of NA1 strains was >1, indicative of evolutionary selection pressure, only two positively selected sites were identified. Both sites are involved in B and T cell epitopes of strain A2 [Norrby et al., 1987; Cane, 1997; Hancock et al., 2003]. These findings are not consistent among different studies, that is, M262 and N273 were also predicted by [Zlateva et al., 2004; Eshaghi et al., 2012], but not confirmed by [Botosso et al., 2009; Tan et al., 2012] using similar methodology. BA-IV strains had higher diversity that appeared to be driven less by positive selection pressure, as the overall dN/dS was 0.45. The two positively selected sites, L219 and T270, have been reported before by [Zlateva et al., 2005; Botosso et al., 2009], while other codons (142, 227, 235, 258, 259, and 311) found to be positively selected in the same studies using stringent criteria were conserved completely in Latvian HRSV-B strains (Fig. 3). There are several possible reasons for inconsistency among the studies. In this study, only one genotype in each group was studied. Other studies have extended over several decades and included other circulating genotypes and extinct lineages. On the other hand, this might also suggest

1979

strain adaptation to varying levels of herd immunity among the different populations. Since the original description in Canada, ON1 strains have been isolated in many parts of the world from children, the elderly and from patients with nosocomial acquisition [Lee et al., 2012]. Prifert et al. [2013] found a strong association of this genotype with intensive care unit admission. However, this association has not been validated in other studies. The ON1-infected patients in this study had lower respiratory tract infections and, even though one of them required oxygen therapy for a brief period of time, they did not require intensive care. Eshaghi et al. [2012] suggested that the 72-nt duplication occurred because of RNA polymerase backtracking at an upstream 7-nt repeat motif and a stem loop in the secondary RNA structure. This is the largest natural insertion in the G gene described to date. Mutations like this are rare and unlikely to appear simultaneously in various locations; therefore, they can be used as a “natural tag” to reconstruct HRSV dissemination [Trento et al., 2010]. To our knowledge, this is the first study that estimates the evolutionary rate, population dynamics and discrete phylogeografic analysis of genotype ON1 based on globally isolated sequences from December 2010 to June 2013. The evolutionary rate determined based on HVR2 sequences (7.92  103 nt/site/ year) was considerably higher than previously reported rates for HRSV-A using the whole genome (6.47  104) and full-length G gene (2.22  103) [Tan et al., 2012]. This rate was also higher than the estimated evolutionary rate of the HVR2 segment of the BA strains (4.7  103) [Trento et al., 2010]. This indicates relaxed selective pressure operating in the duplicated region. The demographic history of ON1 so far has been considerably different than that of the BA genotype, which showed global expansion for the first 10 years after undergoing a 60-nt duplication in the G gene [Trento et al., 2010]. Although the effective population size decreased before the 2012– 2013 season, these data might change when new sequences become available. It is possible that as in the case of the BA genotypes, other mutations besides the duplication are required for optimal adaptation for replication. The mechanism by which duplication provides selective advantage for BA strains is unknown. It has been speculated that it modifies antigenic characteristics of the G protein, allowing escape from neutralizing antibodies. In the case of ON1, this duplication interferes with a region that contains several B- and T-cell epitopes [Cane, 1997] and changes the O- and N-glycosylation properties of the G protein. Discrete phylogeographic analysis suggested spread of genotype ON1 from North America to Western Europe and subsequently to Africa, Asia and Oceania. This hypothesis is supported by several additional pieces of evidence. First, all of ON1 genomes available in GanBank also have E232 ! G and J. Med. Virol. DOI 10.1002/jmv

1980

Balmaks et al.

T253 ! K mutations (for recent global alignment, refer to [Cui et al., 2013]). When NCBI-BLAST was queried for ON1 protein sequences without the duplicated region, none of the retrieved strains had these two mutations; however, E232 ! G was described in an NA1 strain isolated in the same study as genotype ON1 in Canada, suggesting possible ancestry. Second, there was another NA1 strain that was isolated in Canada and in the following season in Italy, Germany, and Latvia (Table II). Third, the calculated tMRCA (August 2010) matches well the estimation by Eshaghi et al. [2012] that this genotype emerged several months prior to winter 2010–2011. To summarize, this study examined the molecular epidemiology of HRSV in Latvia. Overall, that epidemiology reflected global circulation patterns with several local characteristics. The HRSV season started consistently later in Latvia than the average in Europe. High genetic and phenotypic diversity of Latvian strains appeared to be driven only partially by selection pressure. In 2012, a new HRSV-A genotype containing a 72-nt duplication in the G gene emerged in Latvia. By analyzing available international sequences, the global dissemination pattern was reconstructed. Better understanding of HRSV migration could help to develop public health interventions, with the ultimate goal of an HRSV vaccine. ACKNOWLEDGMENTS The authors would like to thank Dr. G. Storch (Washington University School of Medicine, St. Louis, USA) for critical reading of the manuscript. REFERENCES Agrawal AS, Sarkar M, Ghosh S, Chawla-Sarkar M, Chakraborty N, Basak M, Naik TN. 2009. Prevalence of respiratory syncytial virus group B genotype BA-IV strains among children with acute respiratory tract infection in Kolkata, Eastern India. J Clin Virol 45:358–361. Anderson LJ, Dormitzer PR, Nokes DJ, Rappuoli R, Roca A, Graham BS. 2013. Strategic priorities for respiratory syncytial virus (RSV) vaccine development. Vaccine 31:B209–B215. Arnott A, Vong S, Mardy S, Chu S, Naughtin M, Sovann L, Buecher C, Beaute´ J, Rith S, Borand L, Asgari N, Frutos R, Guillard B, Touch S, Deubel V, Buchy P. 2011. A study of the genetic variability of human respiratory syncytial virus (HRSV) in Cambodia reveals the existence of a new HRSV group B genotype. J Clin Microbiol 49:3504–3513. Balmaks R, Kazaks A, Likopa Z, Grope I, Rasnacs O, Gardovska D. 2011. A new tool for molecular monitoring of respiratory syncytial virus in children with lower respiratory tract infections. Riga Stradins University Collection of Scientific Papers 2010:100–105. Balmaks R, Ribakova I, Gardovska D, Kazaks A. 2013. Molecular epidemiology of respiratory syncytial virus during the 2009– 2010 season in Latvia. Arch Virol 158:1089–1092. Bloom-Feshbach K, Alonso WJ, Charu V, Tamerius J, Simonsen L, Miller MA, Viboud C. 2013. Latitudinal variations in seasonal activity of influenza and respiratory syncytial virus (RSV): A global comparative review. PLoS ONE 8:e54445. Botosso VF, Zanotto PMdA, Ueda M, Arruda E, Gilio AE, Vieira SE, Stewien KE, Peret TCT, Jamal LF, Pardini MIdMC, Pinho JRR, Massad E, Sant’Anna OA, Holmes EC, Durigon EL, The VC. 2009. Positive selection results in frequent reversible amino acid replacements in the G protein sene of human respiratory syncytial virus. PLoS Pathog 5:e1000254.

J. Med. Virol. DOI 10.1002/jmv

Cane PA. 1997. Analysis of linear epitopes recognised by the primary human antibody response to a variable region of the attachment (G) protein of respiratory syncytial virus. J Med Virol 51:297–304. Collins PL, Crowe JE. 2007. Respiratory syncytial virus and metapneumovirus. In: Knipe DM, Howley PM, editors. Fields virology. 5th edition. Philadelphia: Lippincott Williams & Wilkins. pp 1602–1646. Collins PL, Melero JA. 2011. Progress in understanding and controlling respiratory syncytial virus: Still crazy after all these years. Virus Res 162:80–99. Connors M, Collins PL, Firestone CY, Murphy BR. 1991. Respiratory syncytial virus (RSV) F, G, M2 (22K), and N proteins each induce resistance to RSV challenge, but resistance induced by M2 and N proteins is relatively short-lived. J Virol 65:1634– 1637. Cui G, Qian Y, Zhu R, Deng J, Zhao L, Sun Y, Wang F. 2013. Emerging human respiratory syncytial virus genotype ON1 found in infants with pneumonia in Beijing, China. Emerg Microbes Infect 2:e22. Dapat I, Shobugawa Y, Sano Y, Saito R, Sasaki A, Suzuki Y, Kumaki A, Zaraket H, Dapat C, Oguma T, Yamaguchi M, Suzuki H. 2010. New genotypes within respiratory syncytial virus group B genotype BA in Niigata, Japan. J Clin Microbiol 48:3423–3427. Delport W, Poon AFY, Frost SDW, Kosakovsky Pond SL. 2010. Datamonkey 2010: A suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics 26:2455–2457. Drummond AJ, Suchard MA, Xie D, Rambaut A. 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol 29:1969–1973. Eshaghi A, Duvvuri VR, Lai R, Nadarajah JT, Li A, Patel SN, Low DE, Gubbay JB. 2012. Genetic variability of human respiratory syncytial virus A strains circulating in Ontario: A novel genotype with a 72 nucleotide G gene duplication. PLoS ONE 7: e32807. Falsey AR, Hennessey PA, Formica MA, Cox C, Walsh EE. 2005. Respiratory syncytial virus infection in elderly and high-risk adults. N Engl J Med 352:1749–1759. Forcic D, Ivancic-Jelecki J, Mlinaric-Galinovic G, Vojnovic G, BabicErceg A, Tabain I. 2012. A study of the genetic variability of human respiratory syncytial virus in Croatia, 2006–2008. J Med Virol 84:1985–1992. Gupta R, Jung E, Brunak S. 2004. Prediction of N-glycosylation sites in human proteins. http://www.cbs.dtu.dk/services/ NetNGlyc/ Hall CB, Weinberg GA, Iwane MK, Blumkin AK, Edwards KM, Staat MA, Auinger P, Griffin MR, Poehling KA, Erdman D, Grijalva CG, Zhu Y, Szilagyi P. 2009. The burden of respiratory syncytial virus infection in young children. N Engl J Med 360:588–598. Hancock GE, Tebbey PW, Scheuer CA, Pryharski KS, Heers KM, LaPierre NA. 2003. Immune responses to the nonglycosylated ectodomain of respiratory syncytial virus attachment glycoprotein mediate pulmonary eosinophilia in inbred strains of mice with different MHC haplotypes. J Med Virol 70:301–308. Hart GW. 1992. Glycosylation. Curr Opin Cell Biol 4:1017–1023. Johnson PR, Spriggs MK, Olmsted RA, Collins PL. 1987. The G glycoprotein of human respiratory syncytial viruses of subgroups A and B: extensive sequence divergence between antigenically related proteins. Proc Natl Acad Sci USA 84:5625– 5629. Julenius K, Molgaard A, Gupta R, Brunak S. 2005. Prediction, conservation analysis, and structural characterization of mammalian mucin-type O-glycosylation sites. Glycobiology 15:153– 164. Khor CS, Sam IC, Hooi PS, Chan YF. 2013. Displacement of predominant respiratory syncytial virus genotypes in Malaysia between 1989 and 2011. Infect Genet Evol 14:357–360. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, Thompson JD, Gibson TJ, Higgins DG. 2007. Clustal W and Clustal X version 2.0. Bioinformatics 23:2947–2948. Lee WJ, Kim Y-j, Kim DW, Lee HS, Lee HY, Kim K. 2012. Complete genome sequence of human respiratory syncytial virus genotype A with a 72-nucleotide duplication in the attachment protein G gene. J Virol 86:13810–13811.

Three-Year Epidemiology of HRSV in Latvia Librado P, Rozas J. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25:1451– 1452. Lozano R, Naghavi M, Foreman K, Lim S, Shibuya K, Aboyans V, Abraham J, Adair T, Aggarwal R, Ahn SY, AlMazroa MA, Alvarado M, Anderson HR, Anderson LM, Andrews KG, Atkinson C, Baddour LM, Barker-Collo S, Bartels DH, Bell ML, Benjamin EJ, Bennett D, Bhalla K, Bikbov B, Abdulhak AB, Birbeck G, Blyth F, Bolliger I, Boufous S, Bucello C, Burch M, Burney P, Carapetis J, Chen H, Chou D, Chugh SS, Coffeng LE, Colan SD, Colquhoun S, Colson KE, Condon J, Connor MD, Cooper LT, Corriere M, Cortinovis M, de Vaccaro KC, Couser W, Cowie BC, Criqui MH, Cross M, Dabhadkar KC, Dahodwala N, De Leo D, Degenhardt L, Delossantos A, Denenberg J, Des Jarlais DC, Dharmaratne SD, Dorsey ER, Driscoll T, Duber H, Ebel B, Erwin PJ, Espindola P, Ezzati M, Feigin V, Flaxman AD, Forouzanfar MH, Fowkes FGR, Franklin R, Fransen M, Freeman MK, Gabriel SE, Gakidou E, Gaspari F, Gillum RF, Gonzalez-Medina D, Halasa YA, Haring D, Harrison JE, Havmoeller R, Hay RJ, Hoen B, Hotez PJ, Hoy D, Jacobsen KH, James SL, Jasrasaria R, Jayaraman S, Johns N, Karthikeyan G, Kassebaum N, Keren A, Khoo J-P, Knowlton LM, Kobusingye O, Koranteng A, Krishnamurthi R, Lipnick M, Lipshultz SE, Ohno SL, Mabweijano J, MacIntyre MF, Mallinger L, March L, Marks GB, Marks R, Matsumori A, Matzopoulos R, Mayosi BM, McAnulty JH, McDermott MM, McGrath J, Memish ZA, Mensah GA, Merriman TR, Michaud C, Miller M, Miller TR, Mock C, Mocumbi AO, Mokdad AA, Moran A, Mulholland K, Nair MN, Naldi L, Narayan KMV, Nasseri K, Norman P, O’Donnell M, Omer SB, Ortblad K, Osborne R, Ozgediz D, Pahari B, Pandian JD, Rivero AP, Padilla RP, Perez-Ruiz F, Perico N, Phillips D, Pierce K, Pope CA III, Porrini E, Pourmalek F, Raju M, Ranganathan D, Rehm JT, Rein DB, Remuzzi G, Rivara FP, Roberts T, De Leo
n FR, Rosenfeld LC, Rushton L, Sacco RL, Salomon JA, Sampson U, Sanman E, Schwebel DC, SeguiGomez M, Shepard DS, Singh D, Singleton J, Sliwa K, Smith E, Steer A, Taylor JA, Thomas B, Tleyjeh IM, Towbin JA, Truelsen T, Undurraga EA, Venketasubramanian N, Vijayakumar L, Vos T, Wagner GR, Wang M, Wang W, Watt K, Weinstock MA, Weintraub R, Wilkinson JD, Woolf AD, Wulf S, Yeh PH, Yip P, Zabetian A, Zheng Z-J, Lopez AD, Murray CJL. 2012. Global and regional mortality from 235 causes of death for 20 age groups in 1990 and 2010: a systematic analysis for the Global Burden of Disease Study 2010. Lancet 380:2095–2128. Martinez I, Dopazo J, Melero JA. 1997. Antigenic structure of the human respiratory syncytial virus G glycoprotein and relevance of hypermutation events for the generation of antigenic variants. J Gen Virol 78:2419–2429. ¨ rvell C, Rafnar B, Norrby E. 1985. Two distinct Mufson MA, O subtypes of human respiratory syncytial virus. J Gen Virol 66:2111–2124. Nikiforova R, Lucenko I, Zamjatina N. 2011. Epidemiology bulletin. Riga: Infectology Center of Latvia. 21 p. Norrby E, Mufson MA, Alexander H, Houghten RA, Lerner RA. 1987. Site-directed serology with synthetic peptides representing the large glycoprotein G of respiratory syncytial virus. Proc Natl Acad Sci USA 84:6572–6576. Ohno A, Suzuki A, Lupisan S, Galang H, Sombrero L, Aniceto R, Okamoto M, Saito M, Fuji N, Otomaru H, Roy CN, Yamamoto D, Tamaki R, Olveda R, Oshitani H. 2013. Genetic characterization of human respiratory syncytial virus detected in hospitalized children in the Philippines from 2008 to 2012. J Clin Virol 57:59–65. Palomo C, Garcı´a-Barreno B, Pe~ nas C, Melero JA. 1991. The G protein of human respiratory syncytial virus: Significance of carbohydrate side-chains and the C-terminal end to its antigenicity. J Gen Virol 72:669–675. Peret TC, Hall CB, Schnabel KC, Golub JA, Anderson LJ. 1998. Circulation patterns of genetically distinct group A and B strains of human respiratory syncytial virus in a community. J Gen Virol 79:2221–2229. Peret TCT, Hall CB, Hammond GW, Piedra PA, Storch GA, Sullender WM, Tsou C, Anderson LJ. 2000. Circulation patterns of group A and B human respiratory syncytial virus genotypes in 5 communities in North America. J Infect Dis 181:1891–1896. Prifert C, Streng A, Krempl C, Liese J, Weissbrich B. 2013. Novel respiratory syncytial virus a genotype, Germany, 2011–2012. Emerg Infect Dis 19:1029–1030.

1981 Reyes M, Eriksson M, Bennet R, Hedlund K-O, Ehrnst A. 1997. Regular pattern of respiratory syncytial virus and rotavirus infections and relation to weather in Stockholm, 1984–1993. Clin Microbiol Infect 3:640–646. Rueda P, Delgado T, Portela A, Melero JA, Garcı´a-Barreno B. 1991. Premature stop codons in the G glycoprotein of human respiratory syncytial viruses resistant to neutralization by monoclonal antibodies. J Virol 65:3374–3378. Shobugawa Y, Saito R, Sano Y, Zaraket H, Suzuki Y, Kumaki A, Dapat I, Oguma T, Yamaguchi M, Suzuki H. 2009. Emerging genotypes of human respiratory syncytial virus subgroup A among patients in Japan. J Clin Microbiol 47:2475–2482. Simoes EAF, Carbonell-Estrany X. 2003. Impact of severe disease caused by respiratory syncytial virus in children living in developed countries. Pediatr Infect Dis J 22:S13–S20. Stensballe L, Devasundaram J, Simoes E. 2003. Respiratory syncytial virus epidemics: The ups and downs of a seasonal virus. Pediatr Infect Dis J 22:S21–S32. Sullender WM, Sun L, Anderson LJ. 1993. Analysis of respiratory syncytial virus genetic variability with amplified cDNAs. J Clin Microbiol 31:1224–1231. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. 2011. MEG A5: Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28:2731–2739. Tan L, Lemey P, Houspie L, Viveen MC, Jansen NJG, van Loon AM, Wiertz E, van Bleek GM, Martin DP, Coenjaerts FE. 2012. Genetic variability among complete human respiratory syncytial virus subgroup A genomes: Bridging molecular evolutionary dynamics and epidemiology. PLoS ONE 7:e51439. Tatochenko V, Uchaikin V, Gorelov A, Gudkov K, Campbell A, Schulz G, Prahl R, Notario G. 2010. Epidemiology of respiratory syncytial virus in children 2 years of age hospitalized with lower respiratory tract infections in the Russian Federation: a prospective, multicenter study. Clin Epidemiol 2:221–227. Trento A, Galiano M, Videla C, Carballal G, Garcı´a-Barreno B, Melero JA, Palomo C. 2003. Major changes in the G protein of human respiratory syncytial virus isolates introduced by a duplication of 60 nucleotides. J Gen Virol 84:3115–3120. Trento A, Viegas M, Galiano M, Videla C, Carballal G, Mistchenko AS, Melero JA. 2006. Natural history of human respiratory syncytial virus inferred from phylogenetic analysis of the attachment (G) glycoprotein with a 60-nucleotide duplication. J Virol 80:975–984. Trento A, Casas I, Caldero
n A, Garcia-Garcia ML, Calvo C, PerezBre~ na P, Melero JA. 2010. Ten years of global evolution of the human respiratory syncytial virus BA genotype with a 60nucleotide duplication in the G protein gene. J Virol 84:7500– 7512. Valley-Omar Z, Muloiwa R, Hu N, Eley B, Hsiao N. 2013. Novel respiratory syncytial virus subtype ON1 among children, Cape Town, South Africa, 2012. Emerg Infect Dis 19:668–670. van Niekerk S, Venter M. 2011. Replacement of previously circulating respiratory syncytial virus subtype B strains with the BA genotype in South Africa. J Virol 85:8789–8797. Venter M, Madhi SA, Tiemessen CT, Schoub BD. 2001. Genetic diversity and molecular epidemiology of respiratory syncytial virus over four consecutive seasons in South Africa: Identification of new subgroup A and B genotypes. J Gen Virol 82:2117– 2124. Wang LF, Collins PL, Fouchier RAM, Kurath G, Lamb RA, Randall RE, Rima BK. 2012. Family paramyxoviridae. In: King AMQ, Lefkowitz E, Adams MJ, Carstens EB, editors. Virus taxonomy: Ninth report of the international committee on taxonomy of viruses. San Diego: Elsevier Academic Press. pp 672–685. White L, Waris M, Cane P, Nokes D, Medley G. 2005. The transmission dynamics of groups A and B human respiratory syncytial virus (hRSV) in England & Wales and Finland: Seasonality and cross-protection. Epidemiol Infect 133:279–289. Yoshida A, Kiyota N, Kobayashi M, Nishimura K, Tsutsui R, Tsukagoshi H, Hirano E, Yamamoto N, Ryo A, Saitoh M, Harada S, Inoue O, Kozawa K, Tanaka R, Noda M, Okabe N, Tashiro M, Mizuta K, Kimura H. 2012. Molecular epidemiology of the attachment glycoprotein (G) gene in respiratory syncytial virus in children with acute respiratory infection in Japan in 2009/2010. J Med Microbiol 61:820–829.

J. Med. Virol. DOI 10.1002/jmv

1982 Zhang Z, Schwartz S, Wagner L, Miller W. 2000. A greedy algorithm for aligning DNA sequences. J Comput Biol 7:203– 214. Zlateva KT, Lemey P, Vandamme A-M, Van Ranst M. 2004. Molecular evolution and circulation patterns of human respiratory syncytial virus subgroup A: Positively selected sites in the attachment G glycoprotein. J Virol 78:4675– 4683.

J. Med. Virol. DOI 10.1002/jmv

Balmaks et al. Zlateva KT, Lemey P, Moe¨s E, Vandamme A-M, Van Ranst M. 2005. Genetic variability and molecular evolution of the human respiratory syncytial virus subgroup B attachment G protein. J Virol 79:9157–9167. Zlateva KT, Vijgen L, Dekeersmaeker N, Naranjo C, Van Ranst M. 2007. Subgroup prevalence and genotype circulation patterns of human respiratory syncytial virus in Belgium during ten successive epidemic seasons. J Clin Microbiol 45:3022–3030.