Ethics Statement

This study was approved by the ethics committee of the Capital Institute of Pediatrics. No written informed consent was obtained because the samples used for this project were collected after routine laboratory testing, and the study was considered less than minimal risk to subjects by the aforementioned committee. The parents or guardians of the children underwent an oral consent process administered by both a clinician and a researcher. Both the clinician and the researcher signed their names in the registration forms to documented this process. The registration forms were obtained. The ethics committee of the Capital Institute of Pediatrics approved this consent procedure.

Sample collection and HRSV antigen detection

Nasopharyngeal aspirates (NPA) were collected from patients who were hospitalized with ARI at the Affiliated Children's Hospital, Capital Institute of Pediatrics during the period from July 2007 to June 2012, and screened for the presence of antigens of 7 common respiratory viruses, including HRSV, influenza A and B, adenovirus and parainfluenza virus I-III, by immunofluorescent assay as a routine test. Indirect fluorescent assay (Chemicon, Temecula, CA, USA) and direct fluorescent assay (Diagnostic Hybrids, Athens, OH, USA) were used for specimens collected from July 2007 to June 2009 and from July 2009 to June 2012, respectively. Briefly, samples were homogenized with 5 ml Hanks solution shortly after collection and followed by centrifuge at 1237×g at room temperature for 5 minutes. The precipitates were used for antigen detection test by following the manufacturers’ instructions and the supernatants were then stored at -80°C for RNA extraction. Some of the supernatants were inoculated into Hep-2 cells for virus isolation before storing at -80°C.

RNA extraction and reverse transcription (RT)

Total RNA was extracted directly from 150 μl of each randomly selected HRSV antigen-positive specimens or supernatant of cell-cultured isolates with Trizol agent (Invitrogen, Carlsbad, CA, USA) by following the manufacturer’s manual, and cDNA was synthesized by using random hexanucleotide primers (Invitrogen, Carlsbad, CA, USA) and moloney murine leukemia virus reverse transcriptase (Invitrogen, Carlsbad, CA, USA).

Multiplex PCR and real time RT-PCR for subgroup-typing, PCR for full length attachment (G) gene

HRSVs were identified as subgroup A or B by using a multiplex PCR with primers GF, GAR and GBR designed according to the sequences of the reference strains in GenBank (Table 1). The PCR for subgroup typing was performed in a 25 μl final volume mixture containing 0.4 μM of each of subgroup A (GAR)- and subgroup B (GBR)-specific reverse primers and 0.4 μM of forward primer (GF), 3 μl of RT products, 0.2 mM deoxynucleoside triphosphates (dNTP), 1 U of DNA polymerase (TransGen Biotech, Beijing, China) and 2.5 μl PCR buffer (TransGen Biotech, Beijing, China) under the following thermocycling conditions: 94°C for 5 minutes, 35 cycles of denaturation at 94°C for 30 seconds, annealing at 49°C for 30 seconds, and extension at 72°C for 30 seconds, followed by a final extension at 72°C for 7 minutes. The amplified products of 269 bp for subgroup A and 353 bp for subgroup B were identified by electrophoresis in 2% agarose gels stained with ethidium bromide and visualized under UV light. Strain Long of prototype HRSV-A and strain CH18537 of prototype HRSV-B were successfully distinguished by this PCR procedure.

Antigen positive while multiplex PCR negative samples were retested by the real time RT-PCR, which primers and probes were designed according to the relative conserved N protein genes for subgroup-typing [33].

The full length attachment (G) genes of both HRSV-A and -B were amplified by PCR with primers S4298 and F164 (Table 1). The PCR was performed in a 25 μl reaction mixture using 0.4 μM of each of forward primer (S4298) and reverse primer (F164), 4 μl RT products, 0.2 mM dNTP, 1.5 mM MgCl₂, 1 U of Platinum^® Taq DNA polymerase (Invitrogen, Carlsbad, CA, USA) and 2.5 μl PCR buffer (200 mM Tris-HCl (pH 8.4), 500 mM KCl) under the following thermocycling conditions: 94°C for 5 minutes, 45 cycles of denaturation at 94°C for 30 seconds, annealing at 49°C for 2 minutes, and extension at 72°C for 30 seconds, followed by a final extension at 72°C for 7 minutes. The products of 1489 bp were analyzed by electrophoresis in 2% agarose gels as aforementioned method. For some samples with weak signal in agarose gel, a 50 μl reaction mixture for PCR was performed in order to meet the quantity requirement for sequencing.

Nucleotide sequencing

The PCR products from full length G genes were sequenced by Invitrogen using cycle sequencing in the forward and reverse direction with primers S4298 and F164.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Monthly distributions of HRSV infection in Beijing during July 2007 to June 2012.

https://doi.org/10.1371/journal.pone.0075020.g001

Sequence analysis

Chromatogram files were inspected using Chromas Lite 2.1 (Thehnelysium, South Brisbane, Australia). Sequence assembly was performed with DNAStar 5.01 (DNASTAR, Wisconsin, USA). Identical sequences were identified with DAMBE software, version 5.3.10 [34].

Phylogenetic analysis

Nucleotide sequences of G genes of HRSV with determined genotypes were obtained from GenBank for reference (Tables S1 and S2). Alignment of sequences from our study with reference sequences was performed on the ClustalW program within MEGA5 software [35]. Nucleotide sequences of HVR2 were aligned on the basis of the aligned predicted protein sequences. The overall mean nucleotide pairwise Jukes-Cantor distances of HRSV-A and -B alignment were 0.05 and 0.06, respectively. Phylogenetic trees of the HVR2 were constructed respectively for HRSV-A and -B alignment by the neighbor-joining method within MEGA5, using nucleotide p-distance for substitution model, complete deletion for gap or missing data treatment and 1000 replicates of bootstrap probabilities for evaluation of confidence estimates. The trees were unrooted.

Glycosylation analysis

Potential N-glycosylation sites were predicted by sequence context of Asn-Xaa-Thr/Ser, where Xaa was any amino acid but proline.

Statistical analysis

Significant differences in positivity rates of HRSV between males and females and between two HRSV seasons were tested with Pearson Chi-Square test using SPSS13.0 software. P-value of <0.05 was considered to be statistically significant.

Nucleotide sequence accession numbers in GenBank

Nucleotide sequences of complete G genes identified in this study have been deposited in the GenBank under accession numbers KC297233 to KC297493 and KC559440.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Unrooted phylogenetic trees of unique sequences of HVR2s of Beijing HRSV-A (A) and HRSV-B (B) strains and reference sequences of identified genotypes.

Trees were constructed using the neighbor-joining method with 1,000 bootstrap replicates by Mega5 software. Only bootstrap values of >70% are shown at the branch nodes. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the p-distance method and are in the units of the number of base differences per site. Beijing strains are colored and indicated by “BJ/”, followed by their isolation numbers. Samples collected from different epidemic seasons are indicated by colored strain names and colored symbols as follows: 2007/2008, ▼; 2008/2009, ◆; 2009/2010, ●; 2010/2011, ■; 2011/2012, ▲. Strain names mentioned in the text are labeled with rectangles to facilitate visualization. New genotype clades are indicated by gray shading The numbers of identical sequences of HVR2 are indicated in parentheses. The genotype assignment is indicated at the right by square brackets. “BA?” refers to HRSV-B strains with a 60 nucleotide duplication that are not assigned to any cluster.

https://doi.org/10.1371/journal.pone.0075020.g002

HRSV distribution

A total of 19,942 NPAs were collected from July 2007 to June 2012, including 7,397 (37.1%) from girls and 12,513 (62.8%) from boys. The gender information for 32 samples was missing. Out of the total NPA samples, 3,160/19,942 (15.8%) were positive for HRSV by antigen detection assay, including 1,079 from girls and 2,074 from boys. The positivity rate of HRSV in boys was significantly higher than that in girls (χ²=13.8, P=0.00).

Samples collected and antigen-positive rate of HRSV in each HRSV season are shown in Table 2. The overall positivity rates of HRSV infection in 2007/2008 and 2011/2012 were significantly higher than each of the other three HRSV seasons. A stronger HRSV activity with early detection and early annual peak of the virus in 2007/2008 season was followed by a weaker 2008/2009 season. This phenomenon has also been reported that a severe and early HRSV activity season is always followed by a weaker late season [36].

The distribution of HRSV in each month is shown in Figure 1. HRSVs were detected from NPAs collected from every month except May and August in 2008, July and August in 2009, July and September in 2010 and July in 2012. HRSV annual epidemics occurred in winter, November to February, in all 5 investigated years. This clear cut seasonality is common in other northern hemisphere countries [25,37-39].

The ages from 16/3,160 HRSV positive patients were unknown. The minimum, lower quartile, median, upper quartile and maximum of ages of HRSV infected patients were 2 days, 2 months, 5 months, 16 months and 194 months, respectively.

HRSV subgroup distribution

Of the 3,160 HRSV antigen-positive NPA specimens, 34 were excluded for further analyses, because these samples were also antigen-positive for other respiratory viruses. A total of 723 samples, including 687 random NPAs and 36 cell-cultured strains isolated from patients who were not overlapped with those 687, were selected for subgroup typing. Among them, 436 and 197 were determined as subgroup A and B by subgroup-specific multiplex PCR, respectively. Twenty-seven and 41 were further determined as A and B by real time PT-PCR, respectively, because they were not able to be typed by multiplex PCR. Twenty-two samples (3.0%, 22/723) were negative for both subgroup-specific multiplex PCR and real time RT-PCR (Table 2). These 22 samples were also negative for full length G gene PCR amplification (data not shown), perhaps because of a non-specific inhibitor of enzyme activity in these clinical specimens. Taken together, 65.9% (462/701) and 34.1% (239/701) of the HRSV were identified as group A and B, respectively. HRSV-A and HRSV-B were co-circulating during these 5 HRSV seasons (Table 2). Subgroup A were predominant during 3 seasons, 2007/2008, 2010/2011 and 2011/2012, while subgroup B were predominant during 2008/2009 to 2009/2010.

Sequence alignments and phylogenetic analysis

Complete sequences of G genes from 190 HRSV-A and 72 HRSV-B were obtained and aligned with representative sequences of identified reference genotypes (Table 2). The phylogenetic trees and strains with identical nucleotide sequences of HVR2 are shown in Figure 2. According to the criteria for a novel genotype proposed by Venter - “sequences clustered together with bootstrap value of 70-100% and with a p-distance of less than 0.07 to all other members in the same phylogenetic cluster” [24], 4 new genotypes, NA3 and NA4 for HRSV-A, and CB1 and BA-C for HRSV-B were identified in this study (Figure 2, gray shading). The bootstrap values of the clusters NA3, NA4, BA-C and CB1 were ≥98%. Furthermore, the maximum p-distances of clusters NA3, NA4, BA-C and CB1 were 0.01, 0.02, 0.02 and 0.04, respectively. As shown in Figure 2, most of the identical HVR2 sequences were from strains isolated in the same HRSV season. However, there were identical sequences isolated from different seasons (e.g. strain BJ/28710 (Figure 2A, green rectangle) and strain BJ/22455 (Figure 2A, yellow rectangle) of genotype NA1, strain BJ/20050 (Figure 2B, magenta rectangle) of genotype BA9). In addition, the HVR2 nucleotide sequence of strain BJ/35320 (Figure 2A, red rectangle), containing a 72 nucleotide duplication, was identical to that of prototype ON1 strain ON67-1210A.

Genotype distribution

Genotype distribution of HRSV in each HRSV season is shown in Table 3. Phylogenetic analysis revealed that 179/190 (94.2%) Beijing HRSV-A strains were classified as NA1 genotype, which was the dominant HRSV-A genotype during the five years. Genotype NA3 was only detected in 2010/2011 season. Genotype NA4 was detected in 2007/2008 and 2011/2012 seasons. One Beijing ON1 strain (BJ/35320) was found from one of the NPAs collected in February, 2012. Most (56/72) of Beijing HRSV-B strains belonged to BA genotype, of which 36 (50.0%) and 10 (13.9%) strains clustered in BA9 and BA-C, while 2 and 14 Beijing HRSV-B strains clustered into SAB4 and new genotype CB1, respectively. The CB1 genotype was detected during four HRSV seasons from 2008/2009 to 2011/2012 (also shown in Figure 2B, gray shading of CB1).

Amino acid analysis

The lengths of all G proteins from Beijing HRSV were deduced based on the complete coding sequences of G genes. The predicted lengths of G proteins from HRSV-B were between 290 and 319 amino acids (Figure 3). It should be noted that a 9 consecutive nucleotide insertion (ACACAAAAA, coding for Thr(T)-Gln(Q)-Lys(K)) occurred next to residue 693 of the G gene (relating to HRSV-B prototype strain CH18537) in the HVR2 in 3 of CB1 strains, BJ/23484, BJ/23183 (Figure 3, black rectangle) and BJ/20265, leading to their G proteins lengthened by 3 amino acids. In addition, although both BJ/F5544 and BJ/20462 (Figure 3, red rectangle) strains clustered in SAB4 genotype, the G protein length of BJ/F5544 was 2 amino acids longer than that of BJ/20462 due to a 6 consecutive nucleotide in-fame deletion in the highly conserved region of G protein of BJ20462 (CCAAAA, residue 475-480 compared to prototype strain CH18537, coding for Pro(P)-Lys(K), data not shown). This deletion was first reported by Zlateva [12] and also seen in genotypes BA7, BA9, BA10 and CB1 of Beijing isolates. The G protein length for HRSV-A was either 297 or 298 amino acids, with the exception of ON1 genotype, which was 321 amino acids due to a 72 nucleotide duplication in HVR2 (Figure 4). However, the G protein length of the Beijing HRSV-A strain BJ/31949 (Figure 4, red rectangle) was 296 amino acids due to a 3 consecutive nucleotide deletion in the highly conserved region of G gene (AAC, residue 571-573 compared to HRSV-A prototype strain A2, data not shown).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Alignments for deduced amino acids of HVR2s from Beijing HRSV-B strains and representative reference strains of identified HRSV-B genotypes.

Amino acid numbering is shown relative to prototype CH18537 (GenBank accession number M17213). Only Beijing strains with unique HVR2 sequences are included in the alignments. Identical residuals, stop codons and gaps are indicated by dots, asterisks and dashes, respectively. Genotype names, followed by the length of G proteins, are shown at the end of the sequences. The G protein length of Beijing strains are deduced from the full length of G protein, while the G protein length of reference strains are given according to prototype of A2. Potential N-glycosylation sites (NXT, where X is not proline) are indicated by gray shading. Different background colors are used to facilitate visualization of different genotypes.

https://doi.org/10.1371/journal.pone.0075020.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Alignments for deduced amino acid of HVR2 from Beijing HRSV-A strains and representative reference strains of identified HRSV-A genotypes.

Amino acid numbering is shown relative to prototype A2 (GenBank accession number M11486). Only Beijing strains with unique HVR2 sequences are included in the alignments. Identical residuals and stop codons are indicated by dots and asterisks, respectively. Genotype names, followed by the length of G proteins, are shown at the end of the sequences. The G protein length of Beijing strains are deduced from the full length of G protein, while the G protein length of reference strains are given according to prototype of A2. Potential N-glycosylation sites (NXT, where X is not proline) are indicated by gray shading. Different background colors are used to facilitate visualization of different genotypes. The 24 amino acid duplication in the genotype ON1 are shown in the rectangle.

https://doi.org/10.1371/journal.pone.0075020.g004

Several genotype-specific substitutions were found in the novel genotypes identified here. Substitutions N255D and S280H were restricted in the HRSV genotype NA3 (Figure 4). The genotype NA4 was characterized by T293S, N242I and P274Q (Figure 4). Genotype BA-C could be defined by M222P, T256A and Q262R. Genotype CB1 could be defined by K224E (Figure 3). Additionally, substitutions R214I and V251M were restricted in strains of CB1 genotype, except for strain BJ/20410 and BJ/28202 (Figure 3, blue and yellow rectangle, respectively). Furthermore, substitution T239I was found only in new genotypes CB1 and BA-C (Figure 3).

Glycosylation analysis

Five major putative N-glycosylation sites, N237, N250, N251, N273 and N294, were identified in HVR2 of HRSV-A (Figure 4, gray shading). Although the first site N237 was conserved in most of genotypes of HRSV-A, it was missing in genotype NA4 and most members of the genotype NA1. However, this site was present in 9 Beijing NA1 strains (Figure 4, blue rectangle). The second site N250 was specific for genotype GA5 and SAA1. All HRSV-A genotypes contained the fifth site N294, except prototype A2 and 6 Beijing NA1 strains (Figure 4, green rectangle). As for HRSV-B (Figure 3), it was interesting to find out that the substitution K258N resulted in an N-glycosylation site in all CB1 strains, which also existed in GB2 genotype and one strain BJ/20240 (Figure 3, green rectangle), but not in any other HRSV-B genotypes. In addition, all BA-C strains lost the N-glycosylation site of N290, which was the same as the prototype CH18537 and genotype BA3.