Epidemics and Frequent Recombination within Species in Outbreaks of Human Enterovirus B-Associated Hand, Foot and Mouth Disease in Shandong China in 2010 and 2011

The epidemiology and molecular characteristics of human enterovirus B (HEV-B) associated with hand, foot and mouth disease (HFMD) outbreaks in China are not well known. In the present study, we tested 201 HEV isolates from 233 clinical specimens from patients with severe HFMD during 2010–2011 in Linyi, Shandong, China. Of the 201 isolates, 189 were fully typed and 18 corresponded to HEV-B species (six serotypes CVA9, CVB1, CVB4, Echo 6, Echo 25 and Echo 30) using sensitive semi-nested polymerase chain reaction analysis of VP1 gene sequences. Phylogenetic analysis based on the VP1 region showed that eight E30SD belonged to a novel sub-genogroup D2; E25SD belonged to a novel sub-genogroup D6; E6SD belonged to sub-lineage C6 and five CVB1SD belonged to subgroup 4C; and B4SD belonged sub-lineage D2. The full viral genomes of the CVB1SD, E6SD, E25SD and E30SD isolates were sequenced. Analysis of phylogenetic and similarity plots indicated that E25SD recombined with E25-HN-2, E30FDJS03 and E4AUS250 at noncontiguous P2A–P3D regions, while E30SD, E30FDJ03, E25-HN-2 and E9 DM had shared sequences in discrete regions of P2 and P3. Both E6SD and B1SD shared sequences with E1-HN, B4/GX/10, B5-HN, and A9-Alberta in contiguous regions of most of P2 and P3. Genetic algorithm recombination detection analysis further confirmed the existence of multiple potential recombination points. In conclusion, analysis of the complete genomes of E25SD, E30SD, CVB1SD and E6SD isolated from HFMD patients revealed that they formed novel subgenogroup. Given the prevalence and recombination of these viruses in outbreaks of HFMD, persistent surveillance of HFMD-associated HEV-B pathogens is required to predict potential emerging viruses and related disease outbreaks.

HFMD usually affects children and has emerged as a significant public health issue in China in recent years. In 2010 and 2011, a total of 3,394,375 cases of HFMD were recorded in Mainland China, accounting for 1,414 deaths, according to the Ministry of Health of the People's Republic of China (http://www.moh.gov. cn/mohjbyfkzj/s3578/201202/54106.shtml). EV-71 and CV-A16 are the major etiological agents of HFMD [12,13,14], but other HEV-A serotypes, CV-A2-8, CV-A9-12, have also been associated with both sporadic infections and outbreaks of HFMD [15,16]. In recent years, some HEV-B serotypes have also been detected in HFMD. A distinct new serotype CVB5 was first reported in 2009 when it was involved in an outbreak of EV-71-associated neurologic HFMD in Linyi City [17]. In addition, other HEV-B serotypes, E1, E6, E11, E13, E25 and E30, were also isolated under the HFMD surveillance network during a 6-month period in 2010, in Henan province, China [18]. E30-and EV71associated HFMD outbreaks were also detected in Guangxi, Southern China in 2011 [19]. However, the epidemiology and molecular characteristics of the HEV-B species associated with HFMD outbreaks in China are still not well known. Moreover, 10-30% isolates from HFMD samples were unable to be characterized in many reports.
In this study, we investigated the presence of enteroviruses in throat swabs from .233 patients in Linyi City in 2010 and 2011, and fully typed 187 isolates from 201 HEVs based on VP1 gene analyses. We present four full-length genomic sequences of the modern HEV-B strains, CVB1SD2011CHN (B1SD), E6SD2011CHN (E6SD), E25SD2010CHN (E25SD) and E30SD2010CHN (E30SD), which represent HEV-B serotypes CVB1, Echo6, Echo25 and Echo30, respectively. We compared the four strains with other related viruses in terms of their partial 59 untranslated regions (UTRs) and P1/P2/P3 regions. Understanding the hereditary and molecular epidemiology of HEV-B will have important implications on epidemiological investigations and the prevention and control of HFMD.

Epidemiological features of enterovirus in Shangdong LinYi
The temporal distribution of the epidemic enterovirus in LinYi City was seasonal, with most cases occurring during the summer from June to September. In 2010-2011, 167 (88.4%) of 189 cases involved individuals younger than 4 years old, and 22 cases (11.6%) involved individuals between 4 and 8 years old. The highest rates of enterovirus-positive samples were from patients younger than 2 years old (75.2% in 2010 and 61.0% in 2011). Of the total 189 isolates, 128 were from males and 61 were from females, giving a male-to-female ratio of approximately 2.15:1. The major clinical symptoms in hospitalized patients with enterovirus-related HFMD were vesicles on the hands, feet, buccal mucosa and tongue, accompanied by fever, vomiting, poor spirit, bilateral sounds of heavy breathing and instability when walking. In addition, at least 50% of EV71-associated HFMD patients also had lung infections, and Cox-As show herpangina symptoms. HFMD patients with HEV-B infection complicated by aseptic meningitis were mostly infected with Echo30 (70%).

Molecular epidemiology of CVB1SD, CVB4SD, E6SD, E25SD and E30SD isolates based on VP1
According to the molecular epidemiological study, genotypes were defined as clusters of related strains with .85% nt sequence identity in the VP1 region [20]. Individual phylogenetic dendrograms for E6, CVB1, CVB4, E25SD and E30SD were drawn based on VP1 sequences from all representative strains responsible for international or domestic epidemics during the past decade. There were four E25 genetic groups, designated A-D (subgroups D1-D6) [11].The E25SD10 sequence detected in the present study belonged to the novel subgenogroup D6 (Figure 2). Four genetic groups were classified for the E30 strains, designated A-D (subgroups D1-D2) [7,21]. The Ukraine and E30 JS2003 strains belonged to subgroup D1, while the E30SD strains isolated in this study and E30Linyi (the latest strains available in GenBank) were categorized as novel subgenogroup D2 with nt sequence identities ranging from 95.6-98.9%, and amino acid identities ranging from 98.9-100% ( Figure 3), which have not been reported previously. The E30 TW, Zhejiang, and GX strains belonged to group C. There were three Echo6 genetic groups designated A-C (subgroups C1-6), and four CVB1 genetic groups designated A-D (subgenogroups D1-6), which had been previously identified [4,9]. Based on this classification, the E6SD sequence detected in the present study belonged to subgroup C6 ( Figure 4) and five B1SD sequences belonged to the sublineage D6 ( Figure 5). Finally, CVB4 was classified into three genetic groups A-C (subgroup C1-2), and two CVB4 sequences formed a cluster in the C2 sublineage ( Figure 6). In addition, phylogenetic dendrograms of the EV71 and CA16 strains indicated that E71 and CA16 formed new subgenogroups C4a and C3c, respectively, as shown in Figures S1 and S2.

Genomic features
We obtained complete genome sequences for the four HEV-B strains, CVB1SD, E6SD, E25SD and E30SD, with lengths of 7,368, 7,395, 7,422 and 7,427 nt, respectively. Their genomic features were analyzed by multiple comparisons with GenBank sequences of other HEV-B enteroviruses that were previously reported as prevalent in China, including 21 prototype strains and 20 modern strains. Pairwise comparisons for CVB1, E6, E25SD and E30SD with the 876-bp VP1 region of other HEV-B strains revealed that the viruses were most similar to prototypes Conn-5B1, Damori E6, E25-JV4, and E30 Bastianni with nt identities of 75.2, 77.2, 81.7 and 82.3%, and amino acid identities of 88.9, 93.6, 92.7 and 95.8%, respectively. VP1 analyses identified the four isolated strains as CVB1, Echo6, Echo 25 and Echo30, respectively. The nt identities of all HEV-B strains in the NSP regions P2 and P3 ranged from 69.1-97%, and the amino acid sequence identities of all HEV-B strains were similar, and ranged from 94.6-98%.

Phylogenetic and recombination analysis of the four HEV-B strains and other HEV-B genomes
The genetic relationships and recombination events between newly isolated HEV-B strains and the other HEV-B strains available in GenBank were investigated by constructing neighborjoining nt phylogenies for the structural protein region (P1 and VP1) (Figure 7a    non-capsid regions indicated that recombination might have occurred during the evolution of these viruses. To address the issue of recombination in more detail, the aligned HEV-B complete genome sequences were further analyzed by examining the similarity among sequences in a sliding window of 400 residues, using the SimPlot program. The sequence of each strain was used as the query sequence to generate 40 separate similarity plots. All HEV-B strain 59NTRs (Figures 8  and 9) were almost equidistant from each other with 85-96% similarity. The VP1 regions of B1SD, E6SD, E25SD and E30SD showed the highest similarities to their prototypes CVB1-conn-5, E6-D9Amori, EV25-JV4, and E30-Bastianni, respectively. CVB5SD, E6SD, E25SD, and E30SD also manifested the highest similarities to their modern strains CVB5-HN, E6-JPN-2011, E25-HN-2, E30FDJ in the 59 half of the genome (ending around nt 4000-5000 within the P2 region). Interestingly, CVB4/JX/10 was highly similar to at least five other strains, E1-JQ979292, CVB5-HN, E6SD, A9-Alberta, and CVB1 (94-98.7%) throughout most of P2 and P3, in which the similarity plot displayed two obvious cross-points at nts 4667 and 5321 among the five strains (Figure 8a and 8b). The similarity plots for E6SD and CVB1 further confirmed that the CVB1SD, CVB5-HN, CVB4/JX/10, A9-Alberta, E1 and E6SD strains were most similar to each other, with high similarity (94-98.7%) extending from nt 4667 or 5321 to the end of 3D (Figure 8c and 8d). In contrast, E25SD and E30SD exhibited a complicated relationship with the other strains E25-HN, FDJS03, E9-DM and E4-AUS250 in the P3 region. E25SD was closely related to FDJS03 in 2C, 3AB, and part of the 3D region (Figure 9a), whereas it shared the highest similarity with E4AUS250 at the 3C-3D junction region. Similarly, E30SD was related to E25-HN-2 in the 3A-3B region, to EV97 at the 3C-3D junction, and to E9 DM in a portion of 3C with over 95% bootstrap values (Figure 9b). In contrast, the regions of similarity between E25-HN2 and FDJS03 differed from E25SD and E30SD according to the temporal order of strain discovery. E25-HN2 was only related to E4AUS250 in 3D with 90% similarity (Figure 9c), whereas E30-FDJS03 was related to E9 in part of 3C with 85% similarity (Figure 9d). The SimPlot results indicated complex mosaic recombination involving CVB1SD11, E6SD11, CVA9-Alberta-2010, COXB5/Henan/2010, CVB4/GX/10, E1-JQ979292, E25SD, and E30SD.

Recombination breakpoints according to GARD
The GARD recombination for putative breakpoints within genomes was confirmed using the HKY85 model and the KH test, which demonstrated significant incongruent topologies before and after each breakpoint (Table 1). Firstly, there was evidence of 10 breakpoints among the alignments of CVB1SD11, B1-1167438 pmmC, E6SD11, E6-FJLY2010, CVA9-Alberta-2010, COXB5/ Henan, CVB4/GX/10, and E1-JQ979292 using AICc; and an improvement of 5,639.97 points was achieved over the model without recombination. Of those, three breakpoints at nts 3487, 4437 and 5109 also indicated significant topological incongruence using the KH test. Ten breakpoints were also detected in the alignments of E25SD, E25-JV4, E30FDJS-3, E21-Farina, E4-AUS-250G, and B5 2000/CSF/KOR; four topology-flanking positions at nts 924, 3430, 5473, and 6393 were significantly discordant by the KH test. Finally, eight breakpoints were detected in E30SD, E30FDJS03-84, E30-Bastianni, EV97 99188/SD/ 1999, and E9DM, three of which (nts 4937, 5549 and 6251) were significantly discordant according to the KH test. The topologyflanking points provide evidence for these significant discordant positions as recombination breakpoints. Overall, the GARD analysis results were consistent with the phylogenies and SimPlot results described above, indicating the occurrence of recombination events in these viruses.

Discussion
Various enteroviral serotypes are spread globally by temporal and regional factors. The dominant causative agent of HFMD outbreaks in China used to be the C4 genotype of EV71 [22]. We previously investigated the variation in and epidemic characteristics of enteroviruses by performing continuous etiological and epidemiological surveillance in Linyi City, Shandong province, China, in 2008 and 2009, and showed that EV71 was the major pathogen of HFMD [13,14]. The results of the present study, however, suggest that the pathogenic spectrum changed in 2010 and 2011, with the most prevalent enterovirus serotype being EV71 in 2010 (81%) but CA16 in 2011 (50%). However, the prevalence of HEV-B species isolated from HFMD cases increased from 3% in 2010 to 17.4% in 2011. In all, the isolation of at least 11 HEV serotypes and 16 HEV-B isolates in the 2011 pathogenic spectrum confirmed that enterovirus species exist not as a single genotype, but as a group of independently evolving genome fragments [23,24]. We tested for the CVB1, CVB4, E6, E25 and E30 serotypes originating from aseptic meningitis outbreaks or other infections in HFMD cases, speculating that HEV-B serotypes may circulate in different clinical diseases in the same region. Our results demonstrated that children younger than 2 years old remained a high-risk population, with an average HEVpositive rate of 68%. A male-to-female ratio of 2.15:1 was found in severe hospitalized HFMD patients, confirming previous investigations. We also found that HEV-B serotypes, especially E30, tended to be associated with severe cases of HFMD complicated by aseptic encephalitis. This represents an increasing etiology in HFMD outbreaks in eastern China and nearby countries, and should thus be a cause of concern in relation to the prevention and control of HFMD epidemics. VP1 clustering or genotyping by phylogenetic analysis can discriminate between lineages within a serotype to identify emerging variants or serotypes, and to investigate the molecular epidemiology of enteroviruses. The E25 phylogenetic tree showed that E25 isolates in China in 1998 (Yunnan and Shangdong provinces) were located in cluster B as a result of temporal evolution. However, nearly all E25 isolates found in China in 2000-2010 were located in cluster D3-D6; the only exception being E25-Mongolia AB239934 isolated in Mongolia in 2003, which appeared in the D3 cluster. The E25SD strain and other E25 strains clustered in sublineage D6 (from Shandong, Henan, and Yunnan provinces) and were most closely related to E25-HN isolated from China in 2008, indicating that the E25-HN isolates may be the ancestors of the 2010 Chinese strains. The sequences of the global E30 isolates can be divided into four genetic lineages (A-D) that display at least 15% diversity between clusters, according to VP1 genotypes [25]. E30 strains from China were first detected in Taiwan in 2001 and were classified in subgroup lineage C1. Phylogenetic tree analyses revealed that all E30 strains To further understand their genetic characters, we sequenced four full-length HEV-B genomes (CVB1SD11, E6SD11, E25SD10 and E30SD10) and aligned them using the ClustalX sequence alignment tool. Based on the updated VP1 sequences in GenBank, the VP1 protein sequences of CVB1SD11, E6SD11, E25SD10 and E30SD10 were 88.9, 93.6, 95.8 and 92.7% identical to those of the corresponding prototype strains, respectively. These are the first complete genomic sequences for recent HEV-B isolates associated with HFMD cases in China. A number of HEV-B serotypes co-circulate thus facilitating spatial  and temporal genetic recombination events [18,28]. Moreover, most circulating HEV-B viruses have incurred multiple recombinations in the NSP region, which commonly results in mosaic enteroviruses [29,30,31,32,33]. The present study showed that recombination occurred during the 2010-2011 HFMD outbreaks. In the structural protein region (P1), E6-SD11 and CVB1-SD11 were most closely related to their respective homotypic prototype strains, whereas in the NSP region, the epidemic E6-SD2011CHN and CVB1-SD2011CHN strains were highly similar to one another, suggesting genetic rearrangement between the E6SD2011 and CVB1SD2011 strains. Interesting, our SimPlot and bootscaning analyses showed that the CVB1SD11, E6SD11, CVA9-Alberta-2010, COXB5/Henan/2010, CVB4/GX/10 and E1-HN serotypes shared sequences throughout most of the P2 and P3 regions, suggesting multiple recombination events. Of those, COXB5/Henan/2010 (HQ998851) and E1-HN (JQ979292) were isolated in 2010 from Henan province, China, and CVB4/GX/10 (JX308222) was isolated in our laboratory in 2010, from Guangxi, China. Their simultaneous presence in neighboring regions would have provided ample opportunity for co-infection and recombination. In present study, we first report the co-occurrence of intertypic recombination within three diverse species, Coxsackie-virus A, Coxsackievirus B and echovirus. E25SD showed high similarity to E25-HN strains in its 59 half, and shared sequences with E30FDJS03 and E4AUS250 in the P2C-P3D region in a discrete manner. Similarly, several recombination sites were also identified in the E30SD genome, in which the regions of similarity were not contiguous among E30SD, E30FDJ, E25-HN, EV97SD and E9 DM strains. Obviously, recombination between E25SD and E30SD only occurred within echovirus species. Thus identical and different HEV-B serotypes co-circulate and recombine within the same geographical area during outbreaks, and recombination produces new variants with modified pathogenic properties. New variants possess serotype-specific capsid protein sequences, but share non-capsid protein sequences of different origins [34]. Given that the analyzed epidemic CVB1SD11 and E6SD11 strains were derived from recombination events, reconstruction of the ancestral state suggests that E1-JQ979292HN-2010 was the donor of the NSP sequences in the genomes of CVB1SD11, E6SD11, CVA9-Alberta-2010, COXB5/Henan/2010, CVB4/GX/10 and E1-JQ979292. Consequently, co-circulation of the CVB1SD, E6SD, E25SD and E30SD strains isolated from HFMD cases with other echoviruses over 2 years resulted in complex mosaic recombination in Shandong and adjacent regions. Attention should therefore In conclusion, the four HEV-B serotypes, E25SD, E30SD, CVB1SD11 and E6SD11 isolated from HFMD patients in Linyi city differed from strains previously isolated in other continents, forming a novel subgenogroup. Given the prevalence and recombination of these viruses in the HFMD outbreaks, persistent surveillance of HFMD-associated HEV-B pathogens is required to allow the prediction of potential emerging viruses and related disease outbreaks.

Ethics Statement
This work received approval from the Clinical Ethics Committee of the Institute of Pathogen Biology, Chinese Academy of Medical Sciences and Peking Union Medical College. Written informed consent was obtained from the next of kin, caretakers, or guardians on the behalf of the minors/children participants involved in our study. The ethics committee specifically approved this consent procedure.
Samples and virus isolation. Throat swabs were obtained from hospitalized HFMD patients at Linyi People's Hospital  (Shandong province, China) during outbreaks in 2010 and 2011. Viral RNA was extracted using an RNeasy Mini kit (Qiagen, Valencia, CA, USA) and stored at 280uC for further use. A total of 233 clinical specimens were subjected to diagnostic reverse transcription-polymerase chain reaction (RT-PCR) generating 350-bp amplicons corresponding to a highly conserved domain in the 59 non-coding (NCR) of HEV. HEV strains were isolated from positive throat swabs by co-culturing in human rhabdomyosarcoma, Vero and buffalo green monkey cells. To avoid viral mixtures, serially diluted samples were prepared and inoculated into cells in 96-well plates, then the most dilute sample that produced a cytopathogenic effect (CPE) was expanded [35]. Viral mixtures were excluded by sequencing the VP1 and 3D genomic regions of individual clones using the TOPOH TA CloningH Kit (Invitrogen, Carlsbad, CA, USA) [36].
Sequence analysis. The amplicons of cDNA sequence fragments were evaluated and assembled manually. The fulllength genome sequences of CVB1SD, E6SD, E25SD and E30SD were aligned using Clustal W (http://www.clustal.org/). All sequences were deposited in the GenBank sequence database (18 HEV-B strains accession numbers: KC411818-KC411825, JQ713883-JQ713888, JX976769-JX976773). The resulting alignments were analyzed using SimPlot software version 3.5.1 (http:// sray.med.som.jhmi.edu/SCRoftware/simplot/) with a sliding window of 400 nucleotides (nt) moving in steps of 50 nt, with each strain as a query in each run. Phylogenetic trees for 3 distinct genomic regions were constructed with the Molecular Evolutionary Genetics Analysis (MEGA) version 5.0 software package (NJ algorithm and Kimura evolution model) [38,39]. Recombination analyses were implemented via the Genetic Algorithms for Recombination Detection [32] program in the DataMonkey software package (http://www.datamonkey.org). Multiple breakpoint detection between the non-recombinant and recombinant models was assessed by comparing the corrected Akaike's Information Criterion (AIC C ) scores. The Kishino-Hasegawa (KH) test was applied to verify whether the adjacent sequence fragments yielded significant topological incongruence. Ancestralstate reconstruction was performed using the neighbor-joining-tree method with DataMonkey software.