Evidence of Recombination and Genetic Diversity in Human Rhinoviruses in Children with Acute Respiratory Infection

Background Human rhinoviruses (HRVs) are a highly prevalent cause of acute respiratory infection in children. They are classified into at least three species, HRV-A, HRV-B and HRV-C, which are characterized by sequencing the 5′ untranslated region (UTR) or the VP4/VP2 region of the genome. Given the increased interest for novel HRV strain identification and their worldwide distribution, we have carried out clinical and molecular diagnosis of HRV strains in a 2-year study of children with acute respiratory infection visiting one district hospital in Shanghai. Methodology/Findings We cloned and sequenced a 924-nt fragment that covered part of the 5′UTR and the VP4/VP2 capsid genes. Sixty-four HRV-infected outpatients were diagnosed amongst 827 children with acute low respiratory tract infection. Two samples were co-infected with HRV-A and HRV-B or HRV-C. By comparative analysis of the VP4/VP2 sequences of the 66 HRVs, we showed a high diversity of strains in HRV-A and HRV-B species, and a prevalence of 51.5% of strains that belonged to the recently identified HRV-C species. When analyzing a fragment of the 5′ UTR, we characterized at least two subspecies of HRV-C: HRV-Cc, which clustered differently from HRV-A and HRV-B, and HRV-Ca, which resulted from previous recombination in this region with sequences related to HRV-A. The full-length sequence of one strain of each HRV-Ca and HRV-Cc subspecies was obtained for comparative analysis. We confirmed the close relationship of their structural proteins but showed apparent additional recombination events in the 2A gene and 3′UTR of the HRV-Ca strain. Double or triple infections with HRV-C and respiratory syncytial virus and/or bocavirus were diagnosed in 33.3% of the HRV-infected patients, but no correlation with severity of clinical outcome was observed. Conclusion Our study showed a high diversity of HRV strains that cause bronchitis and pneumonia in children. A predominance of HRV-C over HRV-A and HRV-B was observed, and two subspecies of HRV-C were identified, the diversity of which seemed to be related to recombination with former HRV-A strains. None of the HRV-C strains appeared to have a higher clinical impact than HRV-A or HRV-B on respiratory compromise.

More than 100 serotypes of HRV are known, which have been classified into two species, HRV-A and HRV-B, according to comparative alignment of nucleotide fragments of VP1 [11,12], VP4/VP2 [13] and 59UTR [14,15], and more recently, on complete genome nucleotide sequences [16]. Moreover, some genomic sequences have been found not to cluster with HRV-A and HRV-B species, which suggests the existence of other species (HRV-C and HRV-D) [16]. A new species of HRV-C was recently identified worldwide by comparative analysis of VP4 or VP4/VP2 genes [7,17,18,19,20,21] and 59UTR [14,15].
However, discrepancies have appeared in the classification of some of the new HRV-A or HRV-C strains, depending on the size and location of the nucleotide sequence in the viral genome and on the phylogenetic methods used for direct analysis of HRV sequences [3,7,14,15,17,18,20,21,22,23].
These recent data underline the lack of knowledge about the biodiversity of HRV strains and their worldwide distribution [14,17]. Moreover, little is known about the characteristics and diversity of HRVs circulating in a given area in a short period of time. In the present study, we looked for HRVs in a 2-year collection of nasopharyngeal swabs (NPSs) of children with ARI visiting a district hospital in Shanghai, and compared sequences in two regions previously defined for genetic classification of HRV serotypes [13,15]. Our study showed a high diversity of HRV species and genotypes, and a prevalence of the novel HRV-C species in NPSs of children with bronchitis and pneumonia. This biodiversity appeared to result partly from recombination events in the 59UTR, between HRV-C strains and those close or similar to HRV-A species, which led to the suggested classification of HRV-C into at least two subspecies.

Identification and typing by phylogenetic analysis of 66 HRVs in NPSs from children with ARI
Eight hundred and twenty-seven samples were collected from a group of children consulting the Shanghai Nanxiang Hospital during a 2-year period, and tested for 17 respiratory viruses using a multiplex RT-PCR (mRT-PCR). Sixty-four samples (7.7%) were positive for HRV, according to the length of the amplified fragment in the VP4/VP2 region visualized on agarose gel (data not shown). A larger fragment of 924 nt, including part of the 59UTR (starting at nt 163) and the VP4/VP2 genes (ending at nt 1086), was amplified and cloned into plasmid vectors for genetic analysis (Table 1). Only one sample, N1, could not be amplified   and was amplified in two steps in the 59UTR (nt 163-552) and in the VP4/VP2 region (nt 528-1086), respectively. Analyses of different clones for each sample allowed the identification of multiple infections: sample N16 was shown to contain one HRV-A and one HRV-B (N16A and N16B, respectively), while N58 contained one HRV-A and one HRV-C (N58A and N58C, respectively) (see below). In order to further characterize the 66 HRVs identified in the 64 samples, the nucleotide sequences located in the 59UTR (285 nt) and the VP4/VP2 genes (420 nt) were chosen to allow comparative alignment with sequences of reference serotypes and field strains available in GenBank.
To characterize and classify further the Chinese HRV strains, 59UTR sequences were considered (Table 3, Fig. 2). They were compared to the 59 UTR of all 101 reference HRVs, to those of 26 new strains identified in children with respiratory illness in Wisconsin (indicated as W) [15], and to those of other recently identified HRV-C strains [14]. Pairwise nucleotide divergence between the three HRV species was 0.7-64.3%, and a limit of ,9% divergence between genotype pairs was chosen for similar genotype assignment in one species [15]. New genotypes were identified when they had 9-30% pairwise nucleotide divergence from the nearest serotype in the same species (Table 3). Fifty-five HRVs shared .94.4% nucleotide identity with strains already identified, and 11 HRVs showed 9.5-20.9% nucleotide divergence with the nearest known HRVs. They may represent newly discovered genotypes. These strains clustered with HRV-A (N6) or HRV-C species (N4, N8, N21, N62, N63, N67, N24, N25, N28 and N32) ( Table 3).
Most surprisingly, 20 of the 34 strains classified as HRV-C by comparative analysis of VP4/P2 sequences ( Table 2) were related more closely to HRV-A strains when their 59UTRs were analyzed, and showed incongruent clustering in phylogenetic trees (Figs. 1 and 2). The nucleotide sequences in the 59UTR of these strains were related closely to those of formerly identified QPM, NAT001, NAT045, C024, C025 and C026 HRV-C strains, and to some of the W strains recently identified as HRV-A [15] (Fig. 2). However, our 20 strains clustered together with these six strains in two major branches in the phylogenetic tree constituted a subspecies of HRV-C called HRV-Ca (Table 3, Figs. 2 and 3). The fourteen other HRV-C strains formed another unique branch of HRV-C subspecies, called HRV-Cc, which clustered differently from other species of HRV-A and HRV-B, and from subspecies HRV-Ca in the phylogenetic tree based on 59UTR sequences (Table 3; Figs. 1 and 2). These strains were related closely to some W strains of HRV that were classified as HRV-C [15] (Fig. 2, Table 3)

Recombination in the 59 UTR
In order to characterize more precisely the differences observed in the 59UTR between HRV-Ca and HRV-Cc subspecies, and to localize possible recombination sites in the 59UTR of the genome of HRV-Ca subspecies, bootscanning and similarity plot analyses were conducted in the gene fragment of 868 nt that included the 59UTR and adjacent capsid genes. HRV-Ca nucleotide sequences were scanned against sequences of N10, R16 and R52, which are considered as representative strains of HRV Cc subspecies and HRV-A and HRV-B species, respectively. Stretches of nucleotide sequences that were closer to HRV-A (R16) than to HRV-Cc (N10), flanked by sequences related to HRV-C could be detected in the 59UTR of HRV-Ca strains, as exemplified with HRV-Ca N25 strain ( Fig. 3a and b). These HRV-A-related nucleotide stretches were thus flanked by putative recombination sites. These sites were located differently among the HRV-Ca strains, delimiting HRV-A-related stretches that ranged from 150 to 400 nt in length (Fig. 3c). While variable among the strains, the identified recombination sites were all located inside the 59 UTR and none of them was identified in the downstream VP4/VP2 coding sequence. HRV-A-related nucleotide sequences and putative recombination sites were also found in the 59UTR of the previously described HRV-C strains C024, C025, C026, NAT001, NAT045 and QPM (Fig. 3). The results corroborated the clustering observed in the phylogenetic tree based on 59UTR sequences (Fig. 2), since strains gathered in the same HRV-Ca subcluster (for example N24, N25, N28 and N32, or N4, N7, N8, N21, N36 and N46) displayed the same recombination pattern. These subclusters revealed different recombinant lineages, each of which originated from independent recombination events.

Comparative analysis of full-length genomes of HRV-A, HRV-B and HRV-C species
In order to further characterize the genome of the HRV-Cc subspecies, for which no full-length sequence was yet available, we sequenced the remaining genes that covered the whole coding sequence and 39UTR of N10 strain, which was chosen as the representative of this subspecies (Tables 2 and 3). The full-length N10 genome sequence was compared to those of the HRV-Ca subspecies strains C024, C025 and C026, and to that of N4 strain, which was sequenced as the representative of the HRV-Ca subspecies. The genome sequences were also compared to those of the HRV-A strains N13 and R44, and to the HRV-B strains R14 and R52 (Table 4). The full-length nucleotide sequence of N10 strain contained 7111 nt, excluding the poly(A) tract, which was shorter than sequences from HRV-A and HRV-B strains, but similar to those of HRV-Ca strain N4 and other related strains (C024, C025 and C026). The 2144 aa lengths of the polyprotein and of each of the individual proteins of N10 were slightly different from those of HRV-A and HRV-B species, but similar to those of other HRV-C strains. The most divergent amino acid length between HRV was observed for the VP1 protein that was shorter in HRV-C species ( Table 4). The unique putative cleavage (M/S) site between VP4 and VP2 protein identified previously for QPM, C024, C025 and C026 strains [20] was also observed for N10 and N4 strains. It was different from those of the HRV-A strains N13 and R44 (Q/S), and from those of the HRV-B strains R14 and R52 (N/S) (data not shown).
Alignment of the VP1 amino acid sequence of HRV-Cc strain N10 with those of other HRV-A and HRV-B species and HRV-Ca subspecies, designated in Table 4, showed structural features typical of HRV-C species [16,20,25] (data not shown). In particular, footprints including deletions in the BC, DE and HI loops and conserved amino acids potentially involved in Inter-Cellular Adhesion Molecule 1 (ICAM-1) receptor binding [7,11,20,24] were conserved within the HRV-C species (data not shown). Bootscanning and similarity plot analysis conducted on the genomic sequences of N4 (HRV-Ca), N10 (HRV-Cc), R16 (HRV-A) and R52 (HRV-B) confirmed that N4 featured a 59UTR sequence that was related to the R16 sequence (stretch I), followed by a capsidic sequence related to the N10 sequence. N4 nonstructural sequence (2A to 39UTR) was related more closely to N10 than to R16 and R52 sequences. However, in stretch II (nt 3300-3500 according to N4 numbering), N4 strain (HRV-Ca) was closer to R16 (HRV-A) than to R52 (HRV-B) or N10 (HRV-Cc), which resulted in high bootstrap values between N4 and R16 2A sequences (Fig. 4). This may have been the result of a recombination event that occurred in the 39 part of the 2A-encoding sequence of the parental strain N4, and which involved an HRV-A strain. Nevertheless, the HRV-A parental strain or ancestral strain could not be identified since the closest HRV-A 2A nucleotide sequence available was ,80% identical to that of N4 in this region.
In contrast from nt 6,550 to the 39 end (stretch III in Fig. 4), the N10 strain genome was found to be closely related to that of N4, with nucleotide identity .98%. This result is corroborated in Figure 5, which shows a phylogenetic analysis of the 39UTR sequences of N10 and N4 compared to those of HRV-Ca subspecies and HRV-A and HRV-B species. This suggests that N4 and N10 strains share a common recent ancestor through recombination.   1 and 2) and on local recombination in 59UTR (see Fig. 3). b Strains closely related with more than 93% identity. c These strains clustered differently when based on VP4/VP2 sequences (see Table 2 and Figs. 1 and 2). doi:10.1371/journal.pone.0006355.t003 Clinical outcome from HRV strains isolated from pediatric outpatients Among the pediatric patients, 46 were males and 18 females, and their age ranged from 5 months to 14 years. The majority of HRV infections were diagnosed between 2 and 6 years of age (84.6%). Bronchitis (73.4%) and pneumonia (26.6%) were highly prevalent in children with comparable incidence in HRV-A and HRV-C infections (Table 5). Moreover, the ratio of pneumonia over bronchitis (36.2%) was comparable to that in the whole cohort of 827 children (40.7%). Only one child among the 64 HRV-positive patients had asthma and was co-infected with HRV-C, influenza A virus (IAV) and respiratory syncytial virus (RSV) ( Table 5), whereas 100 of the 827 patient were diagnosed with asthma. HRVs were isolated throughout the 2 years, with a predominance of HRV-C viruses in the cold season (Table 5). Interestingly, different HRV genotypes were detected within the same period (for example, N1 and N4, N9 and N11/N12, N44/N48 and N51, and N55/N56 and N62/N67), with a larger diversity and distribution of individual or paired HRV-A genotypes compared to HRV-C strains, which clustered in closely-related genotypes (Fig. 2, Tables 2 and 3). Conversely, N4 and N21 strains of samples collected at 10 months interval showed 99.8% identity (Fig. 2).
Single HRV infection was diagnosed in 42 children and coinfections were identified in 22 patients (Table 5), with 17 double and five triple infections. The viruses most often identified in HRV co-infection were RSV (six cases) and human bocavirus (HBoV; four cases), and two patients were co-infected with HRV, HBoV and RSV (Table 5). There was no difference between HRV-Ca or HRV-Cc subspecies and any of the clinical or epidemiological data (data not shown).

Discussion
In this report, we looked for HRVs in a 2-year collection of NPSs from children with ARI visiting a district hospital in Shanghai, and found a high diversity of HRV strains that belonged to different species and genotypes. We characterized by RT-PCR and sequenced 66 HRVs, among them 27 HRV-A, five HRV-B, and 34 HRV-C strains. When sequencing the VP4/VP2 region of the HRV genome, several recent studies have identified new strains of viruses from children and adults with ARI, asthma, or otitis, which are clustered differently from HRV-A and HRV-B, and have been classified into a novel HRV-C species [7,8,17,18,19,20,21,25,26,27,28]. Other groups have also identified novel HRV-C strains by sequencing the VP1 gene [29] or the 59UTR [14,15] . The different sizes and locations of the regions amplified in the HRV genomes renders difficult comparative genetic analysis. Recently, Palmenberg et al. (2009) have finalized the full-length genome sequences of all HRV-A and HRV-B reference strains, and identified structural features of these two species and the novel HRV-C species [16]. In our study, we identified 34 HRVs (51.5%) that clustered differently from HRV-A and HRV-B in a phylogenetic tree that was established on the basis of VP4/VP2 sequences, which were related to recent strains classified in the novel HRV-C species (Fig. 1, Table 2). Fourteen HRV-C strains (41.2%) segregated from the other 20 strains (58.8%) that were closely related to HRV-A in their 59UTR sequence (Fig. 2). This led us to propose a classification of two HRV-C subspecies, HRV-Cc and HRV-Ca. In previous studies  [15]. These strains clustered with our field strains within the HRV-Cc subspecies. Moreover, 17 strains that clustered with HRV-A, and had 12-35% pairwise nucleotide divergence from the nearest reference serotype [15], clustered within the two major branches of HRV-A and HRV-Ca strains (Fig. 3). Therefore, we cannot ensure that some of the 17 strains were HRV-A or HRV-Ca strains. Kiang et al. (2008) have identified five novel HRVs out of 24 clinical samples (20.8%), which segregated from HRV-A and HRV-B, and were classified as HRV-C, and three additional strains (12.5%) that also clustered with QPM, C024, C025, C026, NAT001 and NAT045 [14] (Fig. 2). However, the field HRV strains of these previous studies were sequenced using a 59UTR that did not match fully our sequence and that of Lee et al. (2007) [15], and could not be included in the present study for comparative analysis. Interestingly, the five strains identified in California in 2007 [14] and N42 and N45 from our study were closely related to strain W37 isolated in Wisconsin in the late 1990s [15], and to NAT001 isolated in the winter of 2004 in California [18], which confirms that similar genotypes of HRV-Ca are widespread [17].
The strains of HRV-C species identified in the present study were characterized by analyzing the 59UTR, VP4, and part of VP2 (Fig. 3). This approach showed the advantages of covering only 59NCR, VP4/VP2, VP1 or 3D genome fragments. Analyzing sequences that covered the 59UTR and the downstream VP4/ VP2 capsid region allowed identification of co-infections when several clones were sequenced, and helped to locate the recombination sites in strains of the HRV-Ca subspecies. Thus, this region of the genome may be useful for building a database of the novel strains that are circulating worldwide.
The genome of HEVs is subject to frequent recombination [30,31,32,33,34,35,36], with interspecies exchanges observed in the 59UTR [37]. Palmenberg et al. (2009) have observed intraspecies recombination in three HRV-A, with structural characteristics and phylogenetic evidence that suggests a novel clade D classification [16]. Tapparel et al. (2009) observed phylogenetic incongruities in 59 UTR, VP1 and 3CD sequences of two clinical isolates of HRV-A related to recombination [38]. We observed incongruent clustering of N12, N44 and N48 strains of HRV-A species in the phylogenetic trees based on the 59UTR or the VP4/VP2 regions of their genomes (Figs. 1 and 2, Table 3), which suggests intraspecies recombination in the 59UTR.
We observed one co-infection with HRV-A and HRV-B (N16A and N16B), one with HRV-A and HRV-C (N58A and N58C), and three co-infections of HRV-A and HRV-B with HEVs that may favor recombination events. Previous comparison of genome sequences between 34 HRVs showed only limited recombination events and a pattern of genetic diversity lower than that observed with other picornaviruses [25]. The presence in HRV-C subspecies of sequences that share 90.5-98.6% identity with HRV-A strains (Table 3) suggests that recombination events occurred between HRV-C and HRV-A. Bootscanning of the 59UTR of HRV-C strains also showed different sites and lengths of recombination (Fig. 3c), which suggested that there were several independent events that led to several groups of HRV-Ca genotypes, which formed clusters in the phylogenetic tree (Fig. 2).
Comparative analysis of the full-length nucleotide sequences of two field strains of different HRV-C subspecies (N4 and N10) with those of other HRV species suggested that multiple interspecies recombination events occurred in the 59UTR and in the NS2A protein gene, and that recombination also occurred in the 39UTR between N4 and a strain close to N10. These findings are in agreement with those observed for other HEVs, for which recombination events in the capsid-encoding sequence are very rare, probably because of structural constraints that restrict the functioning of chimeric capsids [31]. This result appeals for the full-length genome sequencing of the major representatives of the HRV-C species, in order to establish a clear understanding of the evolution and classification of the novel virus into subspecies. Comparison of the coding sequences of N10 HRV-Cc with other strains of HRV-Ca subspecies [20], including our field strain N4, showed high similarities in the lengths of the 11 proteins, their cleavage sites, and the structural features of VP1. These characteristics and the absence of growth in cell culture, noted in our laboratory and by others (data not shown), support the classification of the novel strains into a unique HRV-C subspecies.
Our clinical specimens all originated from NPSs from pediatric outpatients. The remarkable outcome of the study is the large diversity of genotypes that has circulated in a relatively small group of people in a district of Shanghai during a 2-year observation. Although some clusters of similar genotypes in a limited period of time were observed, co-circulation of different genotypes and HRV species and subspecies, and co-infections with two HRV species were observed. The prevalence of the novel HRV-C in our specimens (4.1%) differed from previous studies that associated the prevalence of the novel variant with severe disease outcomes, which ranged from influenza-like illness or infection of the low respiratory tract [17,28] to asthma exacerbation, bronchiolitis, and febrile wheeze [7,8,15,18,20,21,29,39]. All our patients showed bronchitis or pneumonia, with no etiological correlation with any of the species or subspecies of HRV. Only one patient co-infected with HRV-C, IAV and RSV was diagnosed with asthma among the HRV-positive patients (1.6%), whereas 100 of the 827 children had asthma (12%). The difference observed with previous studies, 44.6% [8] and 12% [29] asthma in HRV-positive patients, may be related to the criteria for enrolment. Moreover, none of the patients in our study were hospitalized, which makes comparison with hospitalized children difficult [20,24]. Another criterion to consider in the trend to correlate clinical symptoms with HRV infection is the presence of co-infecting pathogens. In our study, four strains of HBoV and six strains of RSV (17.6%) were identified in association with HRV-C (11.7%). HBoV and RSV are common viruses diagnosed in ARI, which are often associated with HRV [40,41], and HBoV was identified in .50% of children co-infected with HRV [20]. Nevertheless, the incidence of HBoV in ARI and in severe outcomes remains elusive [42]. More studies need to be carried out on large numbers of samples from severe and mild diseases, to identify any obvious role of HRV sequence diversity and association with other pathogens in disease severity. Since a large diversity of recombination in HRVs has become obvious, we must be aware of the occurrence of novel HRVs that may become highly virulent.

Materials and Methods
This study was approved by the ethical committee of Shanghai Nanxiang Hospital and written informed consent was obtained from the parents of the children.

Specimens and viruses
Clinical specimens (n = 827) from NPSs were collected from children under 14 years old, who experienced a lower respiratory tract infection, and who were consulting the pediatric department of Shanghai Nanxiang Hospital during the period October 2006 to October 2008.
Samples that showed positive results for HRV were amplified again using specific primers P1-1F and VP4/2R, located in the 59UTR and VP2 gene, respectively (Table 1). One strain of HRV-C (N1) could not be amplified using the P1 and VP2 extreme primers and was amplified using primers in 59UTR and VP4/ VP2, respectively (Table 1). In brief, 2.5 ml of extracted RNA was mixed with 56 buffer and 0.4 mM dNTPs, 0.2 mM of each of the primers, and 1 ml of enzyme mix, and diethylpyrocarbonatetreated ultrapure water was added to a final volume of 25 ml. Amplification programs included reverse transcription at 50uC for 30 min, inactivation at 95uC for 15 min, followed by 40 cycles at 94uC for 30 s, 50uC for 30 s, 72uC for 70 s, and final extension at 72uC for 10 min. The amplified DNA products were detected by ethidium bromide-agarose gel electrophoresis. The lengths of P1-VP2, VP4-VP2 and P1-P3 amplicons were 924, 559 and 390 nt, respectively. DNA products were extracted from agarose gels by using QIAquick Gel Extraction Kit (Qiagen), and were ligated into pMD20-T vector (Takara Biotechnology, Dalian, China), and at least two recombinant plasmids were sequenced in Biosune Sequence Company and Life Biotechnology in Shanghai, China. Sequences of different clones of N16 and N58 showed identities for either HRV-A or HRV-B strains. More plasmids were sequenced for these strains to confirm that the two patients were originally coinfected with two different HRV species.

Complete genome sequencing
Sequences of three complete genomes of HRV were obtained for strains N4 (reference R3061207002 collected on December 7, 2006), N10 (R3070614001 collected on June 14, 2007) and N13 (R3070719007 collected on July 19 2007). Primers used for the amplification of viral genomes were designed after multiple alignments of sequences from the genomes of different HRVs available in GenBank (Table 1). Overlapping amplified DNA products were obtained after PCR of cDNA that was obtained using oligodT and a Transcriptor High Fidelity cDNA Synthesis Kit (Roche, Mannheim, Germany), following the manufacturer's protocols. Briefly, 10.4 ml viral RNA was mixed with 1 ml oligodT and heated at 65uC for 10 min, and then kept on ice for 2 min. After addition of 4 ml 56buffer, 0.5 ml Protector RNase Inhibitor,   2 ml dNTPs, 1 ml DTT, and 1 ml RT enzyme, the reaction was incubated at 50uC for 1 h, inactivated at 85uC for 5 min, and stored at 220uC. Amplification of a 3D region of N4, N10 and N13 HRV strains was carried out by nested-PCR using Takara EXTaq (Takara Biotechnology) and specific primers (Table 1), for 35 cycles of 30 s at 94uC, 30 s at 55uC, and 70 s at 72uC. To amplify VP1 (upstream of 2A) sequences of N4 and N13 strains, nested PCR was carried out using Takara LATaq with GC buffer I(Takara Biotechnology), specific primers (Table 1), and incubation for 35 cycles of 30 s at 94uC, 30 s at 60uC, and 4 min at 68uC. The fragment VP2-3D of N10 HRV strain was obtained by seminested PCR and specific primers VP2 F, and 3D inner and outer reverse primers (Table 1), using Takara LATaq with GC buffer II, for 35 cycles of 30 s at 94uC, 30 s at 60uC, and 6 min at 68uC. The terminal part of the whole genome was obtained by rapid amplification of cDNA ends using 59/39 rapid amplification of the cDNA kit, following the manufacturer's protocol (Roche). To perform 59 terminal RACE, 4 ml of 56cDNA buffer, 2 ml dNTPs, 1.25 ml specific primer 1 (10 mM), 9.2 ml RNA, 1 ml control primer neo1/rev (12.5 mM), 1 ml control RNA, 1 ml RT enzyme, and 0.6 ml RNase inhibitor (Roche) were mixed and incubated for 55uC for 60 min, followed by inactivation at 85uC for 5 min, and stored on ice. The product was purified using the Qiagen PCR Purification Kit and eluted with 30 ml deionized distilled water. A polyA tail was added to the cDNA, by mixing 9.5 ml DNA with 1.25 ml 106 reaction buffer, 1.25 ml (2 mM) dATP, and after incubation at 95uC for 3 min, the reaction was chilled on ice for 2 min. After addition of 0.5 ml terminal transferase, the reaction was incubated at 37uC for 30 min, inactivated at 70uC for 10 min, and kept on ice. Nested PCR was performed by using the Expend High Fidelity PCR kit (Roche). A mixture of 2.5 ml poly-dA-tailed cDNA, 0.5 ml oligodT-anchor primer 37.5 mM, 0.62 ml SP2 primers (10 mM) (Table 1), 0.5 ml control neo2/rev primer (12.5 mM), 0.5 ml dNTP, 0.35 ml enzyme, 2.5 ml 106 buffer, and 18 ml ddH 2 O was incubated for 40 cycles of 30 s at 94uC, 30 s at 60uC, and 30 s at 72uC. To perform 39 terminal RACE, the method was similar to normal two-step RT-PCR using 3D inner and outer F primers (Table 1).
Sequence alignment, phylogenetic analyses and recombination analysis DNA sequences used for P1-P2 gene analysis were based on HRV-16 nt 178-462 and those used for VP4/VP2 gene analysis were based on HRV-16 nt 626-1045. Multiple sequences were aligned using Clustal X [45]. The multiple-sequence alignment was subjected to phylogenetic analyses using programs in the PHYLIP package (v3.6). Bootstrap analysis was performed using SEQBOOT, with a replicate number of 1000. Then, DNADIST and NEIGHBOR were used to obtain distance matrices with the F84 parameter, and a transition/transversion ratio of 4. Consensus trees were computed by CONSENSE, and then re-rooted with RETREE. The final tree was visualized and edited with MEGA version 4 [46].
Recombination analysis was carried out by using Recombination Detection Program v.3.22. Manual bootscanning was performed by using the Juke-Cantor algorithm and the neighbor-joining method [47], with a window size of 200 nt, a step size of 20 nt and 100 replicates. Pairwise identities between sequences were determined with SimPlot software method [48],with a window size of 200 nt and a step size of 20 nt. Pairwise homology matrices were obtained by using CLC Combined Workbench 3.0 software (CLC bio, Aarhus, Denmark).

Nucleotide sequence accession numbers
The original P1-VP2 sequences described in this study were deposited in GenBank under accession nos. GQ223119 to GQ223136. The VP4-VP2 sequences were deposited under the nos GQ223137 to GQ223181, and the P1-P2 sequences under the nos GQ223182 to GQ223226. The full length genomes sequences of N4, N10 and N13 strains were deposited under the nos. GQ223227, GQ223228, GQ223229, respectively.