Molecular epidemiology and clinical characteristics of herpangina children in Beijing, China: a surveillance study

Sample Collection

Between January and December of 2018, pediatric patients with clinical diagnosis of herpangina were monitored by Beijing Tongzhou Center for Disease Control and Prevention (CDC). Children with typically acute onset symptoms, including pharyngeal pain, fever, and herpes or ulcers on the palatine arch, soft palate, uvula and tonsil, but not on the hand, foot or trunk, were defined as herpangina and recruited in the study. The children without the pharyngeal swabs were excluded.

We collected the general demographic date, clinical characteristics and laboratory results for those children with herpangina met inclusive criteria for further analysis.

Ethics statement

This study was approved by Peking University Institutional Review Board Office, with the approval number of IRB00001052-19005. All samples were collected through infectious disease surveillance system which was implemented based on oral informed consent to ensure compliance. Thus oral informed consents were obtained from parents or guardians before recruitment in this study.

Virus extraction and identification

Pharyngeal swabs were collected for nucleotide identification among children with herpangina. RNA/DNA common extraction was performed with QIAamp MiniElute Virus Spin Kit (Cat.no.57704), according to the instructions of the manufacturer. The final extracts were stored in a 1.5 mL Rnase-Free fresh EP tube at −80 °C.

RNA/DNA extracts were amplified by a set of broad-spectrum primers for the 5′ untranslated region (5′UTR) of EVs. Next, specific primers based on VP1 genes were used for the amplification of specimens which were positive for CV-A4, CV-A6, CV-A10 by nested reverse transcription polymerase chain reaction (RT-PCR). All primers and reaction conditions selected in the experiment were designed according to the previous reports and were listed in Table S1. The amplicons products were visible on gels by electrophoresis, and positive products were sent for bidirectional sequencing using Sanger method.

Sequence analysis

The identified nucleotides sequences were assembled utilizing Lasergene’s DNA SeqMan software (version 7.1.0, DNA Star Inc. Madison, WI, USA). Sequence alignment was conducted in nucleotides BLAST board based on NCBI between the obtained VP1 sequences in our study and homologous sequences from GenBank. All sequences identified in this study were submitted to the GenBank database (MN864893–MN864917 for CV-A4, MN864918–MN864953 for CV-A6 and MN864954–MN864977 for CV-A10).

All VP1 sequences of CV-A4, CV-A6 and CV-A10 in the GenBank were retrieved for the construction of the phylogenetic tree. The accessions and details of included strains were acquired from Table S2. The phylogenetic tree was constructed using the maximum-likelihood method with 1,000 bootstrap replications using MEGA (version 7.0.14, developed by Sudhir Kumar, Koichiro Tamura and Masatoshi Nei). Sequence identity matrixes between different branches were calculated by BioEdit (version 7.13, https://itservices.cas.unt.edu/software/bioedit725).

Statistical analysis

Normally distributed continuous variables were described by mean ± standard deviation (SD), and the median and interquartile range (IQR) was used to describe non-normally distributed variables. Categorical variables were described with frequencies and proportions. Comparison between groups was performed by one-way analysis of variance, chi-square test/Fisher exact test or non-parametric test. All statistical steps were implemented by STATA (version 14.0, developed by statacorp LLC, USA).

Demographics characteristics

In 2018, 184 herpangina children were monitored by CDC in Tongzhou routinely, and two outbreaks involved 6 children were reported. Among the 190 herpangina children enrolled in this study in 2018, the median age of was 4.2 (IQR: 2.8–5.5) years old and 99 (52.1%) were boys. The percentages of children aged 0–3 years old, 3–5 years old, and ≥5 years old were 26.3%, 43.7% and 30.0%, respectively (Table 1). The predominant epidemic peak of herpangina children was in the summer season (from May to August) (Fig. 1A).

Enterovirus serotypes

Among the 190 herpangina children, 131 were detected to be positive for enteroviruses with a detection rate of 69.0%. Among the 131 children infected with EVs, eight serotypes were identified in this study, including CV-A6 (39/131, 20.5%), CV-A4 (25/131, 19.1%), CV-A10 (24/131, 18.3%), CV-A16 (5/131, 3.8%), CV-A2 (4/131, 3.1%), EV-A71(2/131, 1.5%), CV-A5 (2/131, 1.5%) and CV-A8 (1/131, 0.8%), while 29 had unidentified serotypes. The predominant epidemic peak of CV-A6 was July and August. The peak of CV-A4 was in July and the peak of CV-A10 was in June, which was described in Fig. 1B.

Clinical manifestation and laboratory results

The clinical manifestations of herpangina children are shown in Tables 1 and 2. Most of the cases (99.5%) were mild, except for one patient with non-EVs infection who had encephalitis and pneumonia. Fever was developed in 95.3% (181/190) of the herpangina children, and 95.3% (181/190) demonstrated herpes. We compared the clinical manifestations of herpangina children infected with three prevalent EVs (CV-A6, CV-A4 and CV-A10) and no significant differences were found between these among groups. As for the laboratory test results of herpangina children, 68.2% (73/107) children showed high white blood cells (>10 ×10⁹/L), 51.4% (55/107) high percentage of neutrophils (>70.0%), and 48.6% (52/107) low percentage of lymphocytes (<20.0%). An elevated serum C-reactive protein (CRP) level was noted among group and the CRP levels in children infected with CV-A10 and CV-A4 were higher than in those infected with CV-A6 adjusted for age and sex (P = 0.025 and P = 0.078, Fig. 2).

Phylogenetic analysis

In laboratory-confirmed cases of enteroviruses, 22 sequences of CV-A4, 36 sequences of CV-A6 and 24 sequences of CV-A10 were amplified by using specific primers in the VP1 region. Based on comparisons with all available VP1 gene sequences downloaded from GenBank, the global phylogenetic analysis of CV-A4, CV-A6 and CV-A10 was performed.

The phylogenetic tree of CV-A4 was constructed with 186 VP1 sequences and 22 amplified Tongzhou strains, and five groups were clustered (Group A-E) in Fig. 3A. CV-A4 strains from Tongzhou showed 84.9%, 83.4%, 87.6%, 91.9% and 94.0% average nucleotide identities with Group A, Group B, Group C, Group D and Group E, and 96.0% homologies intragroup (Fig. 4A). Group A covered strains circulating globally, which indicated that Group A spread and co-evolved across countries. However, almost all the Chinese strains were concentrated in Group E (94.5%), suggesting that this group has become the predominant serotype responsible for the epidemic of CV-A4 in China. Group E was further divided into six additional subgroups (E1–E6) as shown in Fig. 3B. Among them, all Tongzhou CV-A4 sequences were classified in Group E6, which were close to Shandong strains (2011–2014) and Sichuan strains (2014–2015) of China.

The phylogenetic tree of CV-A6 was constructed using 3470 VP1 sequences. All sequences were clustered into six groups (Group A-F). For CV-A6, 83.3%, 89.4%, 90.9%, 90.0%, 92.6% and 95.4% average nucleotide identities were calculated between Tongzhou strains with Group A, Group B, Group C, Group D, Group E and Group F, and the homology of Tongzhou strains intragroup was 96.5% (Fig. 4B). A majority of Tongzhou CV-A6 strains (35/36) detected in this study were distributed in Group F, which was mostly made up of Chinese strains (86.1%), Spain strains (2.75%) and Thailand strains (2.6%). One amplified strain was located in Group E, mainly including the strains from Fukuoka (2015-2016), Yamagata (2015) and Osaka (2015) of Japan (Fig. 5A). It was assumed that this strain might have been imported from Japan. Group F was divided into five subgroups (F1-F5), and the residual 35 strains of Tongzhou were distributed in Group F5 entirely (Fig. 5B). We further segmented Group F5 into five additional subgroups (F5a-F5e). Among them, one strain belonged to Group F5a, which was intimately related to the Yunnan, Shanghai and Jiangxi strains previously reported in 2015–2017 of China. One strain belonged to Group F5d, which was close to the Chinese strains from 2015 to 2017 (Fig. 5B). The remaining 33 strains were entirely located on two branches (F5b and F5c), as shown in Fig. 5C. Subgroup F5b was composed of 195 Chinese strains in 2012–2018, and one Australian strain in 2017, while subgroup F5c mainly consisted of Chinese strains from 2013 to 2017 and Japanese strains in 2017. The circulation of CV-A6 had a distinct spatiotemporal cluster. We observed that the predominant serotypes were Group B with Group C co-circulating before 2008, but after 2008 the dominating ones were Group E and Group F. Almost all the CV-A6 strains circulating contemporarily were in Group F. Furthermore, Group F1 strains were circulated mainly during 2008–2012, but later, Group F5 strains were the primary epidemic (Fig. S1B).

Based on VP1 gene sequences of CV-A10, 801 VP1 sequences and 24 amplified sequences in this study were constructed for the phylogenetic tree. All the strains were classified into five groups (Group A-E). We analyzed the sources of different branches and found that 98.1% of the sequences of Group E were obtained from China (Fig. 6A). CV-A10 strains from Tongzhou revealed 80.5%, 82.6%, 83.3%, 82.3% and 94.7% average nucleotide identities with Group A, Group B, Group C, Group D and Group E, and 95.8% homologies intragroup (Fig. 4C). The phylogenetic tree showed that Group E was divided into six additional subgroups (E1-E6) in Fig. 6B. All the amplified CV-A10 strains were distributed in two groups in Group E5. Ten Tongzhou strains were clustered in a group which was consist of Beijing strains, Henan strains, Yunnan strains and Guangdong strains in 2016–2018, while the remaining 14 Tongzhou stains were largely close to Yunnan strains (2017–2018) in China. Spatiotemporal characteristics were discovered in the subgroups of Group E, which showed that the ascendant cluster was Group E6 between 2009 and 2013. But it was converted to Group E5 after 2013 (Fig. S1C).