key: cord-0048183-eszp9jnw authors: Zhou, Yonghong; Van Tan, Le; Luo, Kaiwei; Liao, Qiaohong; Wang, Lili; Qiu, Qi; Zou, Gang; Liu, Ping; Anh, Nguyen To; Hong, Nguyen Thi Thu; He, Min; Wei, Xiaoman; Yu, Shuanbao; Lam, Tommy Tsan-Yuk; Cui, Jie; van Doorn, H. Rogier; Yu, Hongjie title: Genetic Variation of Multiple Serotypes of Enteroviruses Associated with Hand, Foot and Mouth Disease in Southern China date: 2020-07-28 journal: Virol Sin DOI: 10.1007/s12250-020-00266-7 sha: 05b1e1b43128d87513528d332f7ffe8a55b26cd4 doc_id: 48183 cord_uid: eszp9jnw Enteroviruses (EVs) species A are a major public health issue in the Asia–Pacific region and cause frequent epidemics of hand, foot and mouth disease (HFMD) in China. Mild infections are common in children; however, HFMD can also cause severe illness that affects the central nervous system. To molecularly characterize EVs, a prospective HFMD virological surveillance program was performed in China between 2013 and 2016. Throat swabs, rectal swabs and stool samples were collected from suspected HFMD patients at participating hospitals. EVs were detected using generic real-time and nested reverse transcription-polymerase chain reactions (RT-PCRs). Then, the complete VP1 regions of enterovirus A71 (EV-A71), coxsackievirus A16 (CVA16) and CVA6 were sequenced to analyze amino acid changes and construct a viral molecular phylogeny. Of the 2836 enrolled HFMD patients, 2,517 (89%) were EV positive. The most frequently detected EVs were CVA16 (32.5%, 819), CVA6 (31.2%, 785), and EV-A71 (20.4%, 514). The subgenogroups CVA16_B1b, CVA6_D3a and EV-A71_C4a were predominant in China and recombination was not observed in the VP1 region. Sequence analysis revealed amino acid variations at the 30, 29 and 44 positions in the VP1 region of EV-A71, CVA16 and CVA6 (compared to the respective prototype strains BrCr, G10 and Gdula), respectively. Furthermore, in 21 of 24 (87.5%) identified EV-A71 samples, a known amino acid substitution (D31N) that may enhance neurovirulence was detected. Our study provides insights about the genetic characteristics of common HFMD-associated EVs. However, the emergence and virulence of the described mutations require further investigation. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (10.1007/s12250-020-00266-7) contains supplementary material, which is available to authorized users. Non-polio Enteroviruses (EVs) are among the most prevalent viruses infecting humans. While most infections are asymptomatic, EVs can cause a broad spectrum of disease from mild to severe and even fatal illnesses. Hand, foot and mouth disease (HFMD) is mainly caused by EVs species A and represents an important public health problem in the Asia-Pacific region. China ranks first worldwide in terms of the case number and mortality of HFMD (Yu and Cowling 2019) , and HFMD outbreaks have occurred annually in the past decade. Currently, no specific antiviral drugs are available for HFMD. Despite the availability and partial roll-out of an EV-A71 vaccine, approved for licensure in China in March 2016, the vaccine has no crossprotective effect with other EVs. Thus the burden of HFMD in China remains high . EVs are nonenveloped, single-stranded, positive-sense RNA viruses (https://www.viprbrc.org/brc/home.spg?de corator=picorna_entero) that encode four structural viral proteins (VPs): VP1, VP2, VP3, and VP4. Among these VPs, the VP1 region of EVs is most often used for genotyping (Harvala et al. 2018) . EV-A71, the most common pathogen of severe HFMD, is divided into seven genogroups, A through G; eleven subgenogroups, B0 through B5 and C1 through C5; and two sublineages, C4a and C4b, based on the VP1 region (Iwai et al. 2009; Zhu et al. 2013; Le et al. 2019) . The prototype strain for EV-A71 is BrCr. Coxsackievirus A16 (CVA16) is divided into two genogroups, A and B; two subgenogroups, B1 and B2; and three sublineages, B1a, B1b and B1c, also based on the VP1 region Zhang et al. 2010; Ganorkar et al. 2017) . The prototype strain for CVA16 is G10. Coxsackievirus A6 (CVA6) is divided into four genogroups, A through D; seven subgenogroups, B1 and B2, C1 and C2 and D1 through D3; and two sub-lineages, D3a and D3b (Song et al. 2017) . The prototype strain for CVA6 is Gdula. For both EV-A71 and CVA16, the important neutralizing/antigenic epitopes within the VP1 protein that are associated with viral infectivity have been established. These include the BC-, GH-and EF-loops and numerous residues (EV-A71: residues 31, 43, 98, 104, 106, 145, 164, 242 and 244; 98, [102] [103] 218, 220, 241, 248, 251 and 295) (Huang et al. 2009; Cordey et al. 2012; Lu et al. 2012; Shi et al. 2013; Lyu et al. 2015; Arthur Huang et al. 2017; Song et al. 2017; Xu et al. 2018; Yao et al. 2018; Le et al. 2019) . In contrast, because CVA6 is mostly associated with mild cases and difficult to culture, the corresponding sites and their properties have yet to be determined. In China, previous HFMD outbreaks were mainly caused by EV-A71 and CVA16; however, CVA6 and CVA10 have emerged worldwide as important pathogens in recent years. Currently, CVA6 is the most commonly detected serotype in HFMD patients in China (Song et al. 2017) . To provide a new perspective on viral diversity and identify new potential targets for drug and vaccine research, we conducted phylogenetic and variable site analyses of the VP1 protein of EV-A71, CVA16 and CVA6 viruses isolated from HFMD patients in China. Before the introduction of the EV-A71 vaccine, a prospective virological surveillance program among children hospitalized with HFMD was established between October 2013 and September 2016. The program included three county-level hospitals and three township-level hospitals in Anhua County, Hunan Province, China . Details of the methods and results of the etiological investigations are described elsewhere . Briefly, 2836 patients were enrolled, from whom 2827 throat swabs, 2583 rectal swabs and 778 stool samples were collected. Generic real-time reverse transcriptionpolymerase chain reaction (RT-PCR) was used to detect all EVs, and specific real-time RT-PCR was used to specifically detect EV-A71, CVA16 and CVA6. If a sample tested positive in the generic real-time RT-PCR but negative in the specific real-time RT-PCR, an additional nested RT-PCR was used to amplify a portion of the VP1 or VP4-VP2 regions for viral identification. To select samples for VP1 sequencing, we first calculated the numbers of patients infected with each of the three serotypes (EV-A71, CVA16 and CVA6) by clinical sample type (stool sample, throat swab or rectal swab). These samples were further grouped by year and month of collection to identify specific time points for genotype switches Dolan et al. 2018 ). Finally, we applied a set of selection criteria to these categorized clinical samples to generate a final subsample for VP1 sequencing. The selection criteria were as follows: (1) For throat swabs, when C 10 samples were available in a month, 10 clinical specimens were randomly selected, and when \ 10 samples were available in a month all positive specimens were included; (2) For rectal swabs and stool samples, when [ 2 samples were available in a month, 2 clinical specimens were randomly selected, and when B 2 samples were available in a month, all positive specimens were included; (3) When both fecal and rectal swab samples were available, fecal samples were prioritized, and when fecal samples were not available, rectal swabs were included; and (4) Only one positive specimen was included per patient. VP1 Sequencing of EV-A71, CVA16 and CVA6 Viral RNA was extracted using a QIAamp Viral RNA Mini Kit (QIAGEN, Hilden, Germany) and recovered in elution buffer, according to the manufacturer's instructions. Nested RT-PCR amplification and sequencing of VP1 from the selected clinical samples were conducted using previously described PCR primers (Supplementary Table S1 ) Zhang et al. 2009; Gaunt et al. 2015) . In short, the first round amplification was performed in a 50 lL reaction volume using a SuperScript III One-Step RT-PCR Kit (Invitrogen, Carlsbad, CA, USA), with annealing temperatures of 50°C for EV-A71, 45°C for CVA16 and 55°C for CVA6 over 40 cycles. CVA16 and CVA6 were amplified in a second round using DreamTaq Green PCR Master Mix (Invitrogen) under the same PCR conditions as the first-round in a 25 lL volume. The expected size of the PCR products was approximately 900 bp. The obtained PCR products were sequenced bidirectionally using a Big Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems Inc., Foster City, CA, USA), using corresponding forward and reverse PCR primers in a 96-capillary 3730xl DNA Analyzer (Applied Biosystems). Sequences were confirmed by BLAST search (https://blast. ncbi.nlm.nih.gov/Blast.cgi). For clarity, we first removed all sequences that were \ 1% divergent from each other at the nucleotide level using the CD-HIT program (http://weizhong-lab.ucsd.edu/cdhit_ suite/cgi-bin/) (Huang et al. 2010) . Recombination was detected using a combination of methods (RDP, GENE-CONV, MaxChi, Chimaera, 3Seq, Bootscan and Siscan) within Recombination Detection Program version 4 (RDP4) with the default settings (Martin et al. 2015) . Recombination was considered present when more than three of the above methods showed significant values (P \ 0.05), and was reconfirmed by phylogenetic analysis. Sequence alignment was carried out using Geneious alignment in Geneious v7.1.3 (https://www.geneious.com). The maximum likelihood (ML) method available in IQ-TREE v1.4.3 was employed to reconstruct phylogenetic trees of the corresponding EVs serotypes (Nguyen et al. 2015) . A bootstrap procedure with 10,000 replicates was used to assess the reliability of the obtained phylogenetic trees. Tree files were visualized using FigTree v1.4.3. VP1 amino acid sequence analysis was conducted using PredictProtein (https://predictprotein.org/) (Yachdav et al. 2014) . Selection pressure was measured by identifing instances of episodic positive selection, at the level of an individual site, using a mixed effects model of evolution (MEME) as previously described (http://www.datamonkey. org) (Murrell et al. 2012) . To identify VP1-region amino acid changes in EV-A71, CVA16 and CVA6, sequences were compared to their respective reference strains (genotype A). All VP1 amino acid sequences obtained in this study were aligned. The sequences obtained in this study were submitted to the NCBI under accession numbers MK306732 through MK307505. In total, 2836 HFMD patients were enrolled in our previous study. Using PCR methodologies, EVs that were eligible for genetic characterization in this study were detected in 2517 (89%) patients ). Among the 18 EVs serotypes identified, CVA16, CVA6 and EV-A71 were the most frequently detected, accounting for 32.5% (n = 819/2517), 31.2% (n = 785/2517) and 20.4% (n = 514/2517), respectively (Fig. 1A) . Additionally, the serotypes were identified as follows: 5.9% (n = 149/2517) CVA10, 3.0% (n = 55/2517) CVA4, 1.6% (n = 40/2517) CVA8, 1.2% (n = 31/2517) CVA2, 0.5% (n = 12/2517) E18, 0.3% (n = 8/2517) CVA5, 0.3% (n = 8/2517) CVB5, 0.2% (n = 6/2517) CVB2, and 0.2% (5/2517) E9. Other EVs serotypes were detected in 2 or fewer patients. According to the selection criteria, we selected 17% of the CVA16 isolates (n = 136/819), 23% of the CVA6 isolates (n = 181/785) and 28% of the EV-A71 isolates (n = 144/ 514) for complete VP1 sequencing. Then, we successfully obtained complete VP1 sequences from 76% of the CVA16 isolates (n = 103/136), 97% of the CVA6 isolates (n = 176/181) and 50% of the EV-A71 isolates (n = 72/ 144). After the removal of identical sequences (where variation was \ 1%), 24 EV-A71, 43 CVA16 and 65 CVA6 VP1 sequences remained for analysis ( Fig. 1B-D , Table 1 ). All sequences fell within 3 subgenogroups: C4a for EV-A71, B1b for CVA16 and D3a for CVA6 (Figs. 2, 3 and 4) . Partial VP1 region sequences of CVA10 (n = 119/149), CVA4 (n = 55/75), CVA8 (n = 38/40), CVA2 (n = 16/31), E18 (n = 10/12), CVB5 (n = 8/8), CVA5 (n = 7/8) and CVB2 (n = 6/6) were obtained (Table 1) . After removing sequences with \ 1% variation, partial VP1 sequences of CVA10 (n = 66/119), CVA4 (n = 29/55), CVA8 (n = 28/ 38), CVA2 (n = 11/16), E18 (n = 8/10), CVB5 (n = 8/8), CVA5 (n = 6/7) and CVB2 (n = 5/6) were included for analysis (Table 1) . Following the calculation of similarity in terms of both nucleotides and amino acids, the results showed little variation (Supplementary Table S2 ). Thus, based on the phylogenetic tree and similarity analysis, these isolates clustered within their 8 subgenogroups: D, I-A, II, D, C, D2, E and 5.3 (Supplementary Figure S1-S8 ). The remaining sequences were compared to reference prototype strains to identify variations at the amino acid level. For the EV-A71 sequences, when compared to the prototype strain BrCr, 30 polymorphic amino acid sites were found. For the CVA16 sequences, when compared to the prototype strain G10, 29 polymorphic amino acid sites were found. For the CVA6, when compared to the prototype strain Gdula, 44 polymorphic amino acid sites were found (Figs. 5, 6 and 7, Supplementary Table S3) . Additionally, we identified amino acid residue variations associated with major neutralization/antigenic epitopes. For the EV-A71 sequences, 87.5% (n = 21/24) of samples exhibited a substitution of aspartic acid (D) to asparagine (N) at residue 31 (D31N). Other substitutions were also found in the EV-A71 sequences: N104D/H (n = 2/24, 8%), R145E (n = 24/24, 100%), G239E (1/24, 4%), E244K (n = 24/24, 100%), S283T (15/24, 62.5%) and A293S/G (16/24, 66.7%) ( Table 2) . For the CVA16 sequences, the substitutions L23I/M/V and I235V were found in 88.4% (n = 38) and 9.3% (n = 4) of the investigated samples (n = 43), respectively. Other substitutions were found in only one CVA16 sequence: P27S, L146X, V147S and T289A (Table 3) . 1.5%), G160N/S (n = 65/65, 100%), T205I (n = 1/65, 1.5%) and Q216H (n = 1/65, 1.5%) ( Table 4 ). Despite the variations described above, an analysis of selection pressure using MEME revealed that none of the protein-coding sequences showed adaptive evolution for EV-A71, CVA16 or CVA6. The present study describes the genetic characteristics of EVs serotypes associated with HFMD in Anhua County, Hunan Province, China, between 2013 and 2016. The VP1 protein plays a major role in the infectivity, replication and virulence of EVs and is the primary target for neutralizing antibodies. As such, monitoring surveillance of changes within the VP1 region of EVs may allow for the early identification of epidemiological and clinical changes, as well as the identification of potential novel vaccine targets. All genogroups (or subgenogroups) of the EVs investigated in this study were consistent with the circulating strains reported in mainland China in recent years (Zhang et al. 2010; Baek et al. 2011; Hu et al. 2011; He et al. 2013; Zhang et al. 2013; Chen et al. 2016; Yang et al. 2016; Song et al. 2017; Zhang et al. 2017; Chen et al. 2018) . We found that despite the complex epidemic patterns of HFMD recorded in Anhua County , each serotype investigated demonstrated long-term, population-level, persistent circulation of specific subgenogroups. Since its first detection in 1998, EV-A71 subgenogroup C4 has persistently circulated in mainland China, while EV-A71 subgenogroup B5 has been reported only sporadically . Notably, the C4a subgenogroup is now endemic to China and has been responsible for the ongoing nationwide HFMD outbreak since 2008. Further, our study indicates that C4a remains the single circulating EV-A71 subgenogroup. This finding is in contrast to the EVs circulation patterns reportrd in other areas, such as Taiwan (China), Japan, Singapore, Thailand and Vietnam, in which EV-A71 subgenogroup replacement has been frequently reported Linsuwanon et al. 2014; Le et al. 2019) . Similarly, the cocirculation of multiple lineages (e.g. B1a, B1b and/or B1c) of CVA16 has been observed in Taiwan (China), Malaysia, Thailand, Australia, Vietnam, and Japan (Zhang et al. 2010; Zhao et al. 2016; Sun et al. 2017) , while subgenogroups have remained stable in mainland China. A previous analysis of CVA16 epidemic strains found that subgenotype B1 has been the dominant lineage in China since 2000 (Sun et al. 2017) . Here, our phylogenetic analysis indicated that, the main circulating CVA16 subgenogroup between 2013 and 2016 was B1b. In contrast to that of EV-A71 and CVA16, the molecular epidemiology of CVA6 detected in this study was similar to that in the broader Asian region. The emergence of CVA6 as a major cause of atypical HFMD worldwide has been associated specifically with genotype D (Song et al. 2017) . HFMD-associated CVA6 subgenogroup D3 was first described following an outbreak in Finland in 2008, followed by other countries across Europe (Bian et al. Virologica Sinica 2015; Song et al. 2017 ). Subsequently, the D3 subgenogroup emerged in Asia, and was described from samples taken during outbreaks in Japan, Vietnam and other countries (Anh et al. 2018) . In this study, we found that the D3a lineage was the predominant virus in China between 2013 and 2016, and no indication of other emerging subgenotypes was observed. Although the identified EV-A71, CVA16 and CVA6 sequences belonged to a single subgenogroup across our study period, novel amino acid changes were also detected when compared to the sequences reported in previous molecular epidemiological studies. To our knowledge, amino acid substitutions in the EV-A71 VP1 region at sites Q48R, A54T/V, N104D/H, G239E, T295I/V/S and T296N have been described only in this study. Furthermore, we also found that the detected frequencies of polymorphic sites at positions S283T and A293S/G of the VP1 region were much higher than those in previous reports (van der Sanden et al. 2010; Le et al. 2019) . Interestingly, the amino acid variation D31N in the VP1 region of EV-A71 was previously reported to be associated with higher virulence and neurotropism (Zhang et al. 2014; Le et al. 2019) . However, while the D31N mutation was present in most EV-A71 samples (21/24, 87.5%) in our study, all patients with EV-A71 were clinically mild cases, which indicates that D31N may not be the only factor or even a major factor determining the virulence of EV-A71. However, it is worth noting that although all enrolled patients in this study were hospitalized, only one enrolled patient had severe HFMD. This is due to the healthcare system structure in China, where rural healthcare insurance provides a larger proportional reimbursement for costs for inpatient care compared to outpatient care. As such, parents of sick children request hospital admission regardless of disease severity and physicians support these requests due to concern over the development of severe HFMD illness in seemingly mild cases. Unfortunately, this process restricted our ability to compare genetic variation between different groups of disease severity in this study. Among CVA16 samples, we also detected novel variations at sites P27S, L146X and V147S of the VP1 region. Further, the previously reported variation at site L23I/M/V was detected at higher frequencies in our samples than in previous reports (Sun et al. 2017; Xu et al. 2018) , which may be associated with local circulation and environmental factors. However, despite these findings, we did not detect a large degree of genetic diversity in CVA16 viruses, which is consistent with previous findings (Zhao et al. 2016) . CVA6 is a newly emerging virus that is often associated with milder cases and outbreaks of HFMD. Due to the emerging status of this virus, the variation analysisof CVA6 was based on known amino acid sequences and studies of EV-A71 and CVA16. Here, we detected ten variations in the VP1 region of 102, 106, 151, 160, 165, 174 and 216) , that are potentially associated with major neutralization/antigenic epitopes. Additionally, CVA6 samples from our study demonstrated higher rates of the amino acid variation G160N/S in the VP1 region than previously reported. Based on the same analysis method employed in this study, Kanbayashi et al. (2017) showed that only 68% (n = 17) of samples contained the G160N/S variation in the VP1 region. In contrast, 100% (n = 65) of our CVA6-positive samples contained this amino acid change. Furthermore, we detected previously unreported variations at sites 96-97, 141, 151, 165, 205 and 216 , which were the main drivers of diversity within the VP1 regions of the CVA6 viruses found in our study. The above variations in the VP1-region amino acids of EV-A71, CVA16 and CVA6 may contribute to future changes in virulence, antigenic properties or genotype switches in circulating EVs in China. However, because only mild clinical cases of HFMD were enrolled in our study, we could not compare genotype variations between HFMD severity groups. Furthermore, this study only investigated recombination in the VP1 region, while other groups of investigators studied the entire viral genome and found indicators of recombination more frequently occurring within nonstructural proteins (e.g., P2 and P3) (Oberste et al. 2004; McWilliam Leitch et al. 2012) . Nevertheless, the data presented here may still contribute to drug target research and multivalent vaccine development. In conclusion, this study has enriched the data on the genetic characteristics of multiple HFMD-associated EVs. EV-A71 subgenogroup C4a, CVA16 subgenogroup B1b and CVA6 subgenogroup D3a were the predominant EVs lineages in Anhua County, Hunan Province, China. This phylogenetic analysis of the main serotypes of EVs causing HFMD adds to our knowledge about enteroviral evolution and circulation. However, further research regarding amino acid variations and their effects on virulence, antigenic shifts and genotype switches of EVs is required. 18XD1400300), the Li Ka Shing Oxford Global Health Programme (No. B9RST00-B900.57), the Chinese Preventive Medicine Association (No: 20101801) . JC is supported by CAS Pioneer Hundred Talents Program. Author Contributions HJY conceptualized, designed and supervised the study. KWL, QHL, PL, MH and SBY coordinated and participated in data collection. YHZ, LLW and QQ performed the experiments. GZ, NTA, NTTH, XW, TTL and JC provided laboratory technical support. YHZ and LVT wrote the first draft of the manuscript. HJY, HRD and LVT revised the manuscript. All authors contributed to review and revision and approved the final manuscript as submitted and agree to be accountable for all aspects of the work. Conflict of interest HJY has received investigator-initiated research funding from Sanofi Pasteur, GlaxoSmithKline, and Yichang HEC Changjiang Pharmaceutical Company, none of which is related to HFMD and enteroviruses. Animal and Human Rights Statement This study was approved by the ethical review committees at the Chinese Center for Disease Control and Prevention, World Health Organization Regional Office for the Western Pacific, and the School of Public Health, Fudan University . Informed verbal consent was obtained from all individual patients' parents/guardians included in the study . Emerging coxsackievirus a6 causing hand, foot and mouth disease Epitope-associated and specificity-focused features of ev71-neutralizing antibody repertoires from plasmablasts of infected children Epidemics of enterovirus infection in Chungnam Korea Coxsackievirus A6: a new emerging pathogen causing hand, foot and mouth disease outbreaks worldwide Comparative genetic analysis of VP4, VP1 and 3D gene regions of enterovirus 71 and coxsackievirus A16 circulating in malaysia between Genomic characteristics of coxsackievirus A8 strains associated with hand, foot, and mouth disease and herpangina Multiple transmission chains of coxsackievirus A4 cocirculating in china and neighboring countries in recent years: phylogenetic and spatiotemporal analyses based on virological surveillance Monitoring antigenic variations of enterovirus 71: implications for virus surveillance and vaccine development Identification of site-specific adaptations conferring increased neural cell tropism during human enterovirus 71 infection Mechanisms and concepts in RNA virus population dynamics and evolution Genetic characterization of enterovirus strains identified in hand, foot and mouth disease (HFMD): emergence of B1c, C1 subgenotypes, E2 sublineage of CVA16, EV71 and CVA6 strains in india Spectrum of enterovirus serotypes causing uncomplicated hand, foot, and mouth disease and enteroviral diagnostic yield of different clinical samples Genetic characterization of human coxsackievirus A6 variants associated with atypical hand, foot and mouth disease: a potential role of recombination in emergence and pathogenicity Recommendations for enterovirus diagnostics and characterisation within and beyond europe Emergence, circulation, and spatiotemporal phylogenetic analysis of coxsackievirus A6-and coxsackievirus A10-associated hand, foot, and mouth disease infections from Complete genome analysis of coxsackievirus A2, A4, A5, and A10 strains isolated from hand, foot, and mouth disease patients in china revealing frequent recombination of human enterovirus a Reemergence of enterovirus 71 in 2008 in Taiwan: dynamics of genetic and antigenic evolution from 1998 to CD-HIT suite: a web server for clustering and comparing biological sequences Genetic changes of coxsackievirus a16 and enterovirus 71 isolated from hand, foot, and mouth disease patients in Toyama Impact of coxsackievirus a6 emergence on hand, foot, and mouth disease epidemic in osaka city Molecular epidemiology analysis of enterovirus 71 strains isolated in dak lak Emerging enteroviruses causing hand, foot and mouth disease Epidemiology and seroepidemiology of human enterovirus 71 among Thai populations An atypical winter outbreak of hand, foot, and mouth disease associated with human enterovirus 71 Enterovirus 71 disrupts interferon signaling by reducing the level of interferon receptor 1 Crystal structures of yeastproduced enterovirus 71 and enterovirus 71/coxsackievirus a16 chimeric virus-like particles provide the structural basis for novel vaccine design against hand-foot-and-mouth disease RDP4: detection and analysis of recombination patterns in virus genomes The association of recombination events in the founding and emergence of subgenogroup evolutionary lineages of human enterovirus 71 Detecting individual sites subject to episodic diversifying selection IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies Evidence for frequent recombination within species human enterovirus b based on complete genomic sequences of all thirty-seven serotypes Evaluation of different clinical sample types in diagnosis of human enterovirus 71-associated hand-foot-and-mouth disease Molecular phylogeny of modern coxsackievirus A16 Identification of conserved neutralizing linear epitopes within the VP1 protein of coxsackievirus A16 Persistent circulation of coxsackievirus A6 of genotype D3 in mainland of china between Epidemiological characterizations, pathogen spectrum and molecular characteristics of coxsackievirus A16 from patients with HFMD in Yantai Evolutionary trajectory of the VP1 gene of human enterovirus 71 genogroup b and c viruses Genetic characteristics of the P1 coding region of coxsackievirus A16 associated with hand, foot, and mouth disease in China Predictprotein-an open resource for online prediction of protein structural and functional features Two genotypes of coxsackievirus A2 associated with hand, foot, and mouth disease circulating in china since Poorly neutralizing polyclonal antibody in vitro against coxsackievirus a16 circulating strains can prevent a lethal challenge in vivo Remaining challenges for prevention and control of hand, foot, and mouth disease An outbreak of hand, foot, and mouth disease associated with subgenotype c4 of human enterovirus 71 in shandong Molecular evidence of persistent epidemic and evolution of subgenotype B1 coxsackievirus A16-associated hand, foot, and mouth disease in china Complete genome analysis of the C4 subgenotype strains of enterovirus 71: predominant recombination C4 viruses persistently circulating in china for 14 years The variations of VP1 protein might be associated with nervous system symptoms caused by enterovirus 71 infection Molecular characterization of a new human coxsackievirus B2 associated with severe hand-foot-mouth disease in yunnan province of china in 2012 Characterization of VP1 sequence of coxsackievirus A16 isolates by bayesian evolutionary method Phylogenetic analysis of enterovirus 71 circulating in beijing Acknowledgements We thank staff members of the Anhua County-, Yiyang Prefecture-, and Hunan Provincial-level departments of health for providing assistance with administration and data collection; staff members at the Anhua County-, Yiyang Prefecture-, and Hunan Provincial-level CDCs and six study hospitals (Anhua People's Hospital, Anhua Second People's Hospital, Anhua Hospital of TCM, Tianzhuang Township Hospital, Jiangnan Township Hospital, and Qingtang Township Hospital) for providing assistance with field investigation, administration and data collection.