Introduction

Hand, foot and mouth disease (HFMD), which is mainly caused by some members of the species Enterovirus A, is a common childhood infectious disease [1, 2]. Enterovirus A is one of 15 species of the genus Enterovirus in the family Picornaviridae (http://www.picornaviridae.com/enterovirus/enterovirus.htm). This species includes viruses of 25 types, of which coxsackievirus A2 (CVA2), CVA3, CVA4, CVA5, CVA6, CVA8, CVA10, CVA12, CVA14, and CVA16 and enterovirus A71 (EV-A71) are associated with a variety of human diseases, such as HFMD, herpangina, encephalitis, paralysis, myelitis, meningitis, upper and lower respiratory diseases, and gastroenteritis [2].

HFMD has been one of the major public health problems in mainland China since 2008 when HFMD was listed as a notifiable disease. According to the national notifiable disease reporting system, a total of 11,501,230 HFMD cases, including 955 lethal cases, were reported to the Chinese Center for Disease Control and Prevention (China CDC) in mainland China during 2014-2018, with approximately 2.3 million reported cases each year. The data indicate that although most HFMD patients have only mild symptoms, some develop severe complications that can be fatal. We previously reported that EV-A71 was the main pathogen of severe HFMD in Shenzhen, China, from 2009 to 2013 [3]. Other EV-A types, including CVA2, CVA4, CVA5, CVA6, CVA8, CVA10, CVA14 and CVA16, can also cause severe HFMD [3,4,5,6,7,8,9,10]. Some members of the species Enterovirus B, such as CVA9, CVB1-5, echovirus 3 (E3), E6, E7 E9, E11, E14, E16, E18, E20, E30 and E33, are also frequently associated with severe complications, especially encephalitis and aseptic meningitis [11,12,13,14,15,16,17].

Numerous studies have focused on the epidemiology, functional genomics, and pathogenesis of EV-A7 and vaccine development [18,19,20,21]. It is known that multiple sites in the genome of EV-A71 are associated with its virulence or clinical phenotype [19, 20]. Although a commercial EV-A71 vaccine was available in mainland China in 2016, little is known about the molecular epidemiology of the other enteroviruses (EVs) related to severe HFMD, especially the relationship between genetic features of the virus and disease severity. The aim of this study was to investigate the epidemiology and molecular characteristics of EVs associated with severe HFMD in Shenzhen, China, from 2014 to 2018 before and after the introduction of the EV-A71 vaccine. We focus on severe complications in patients infected with different EVs, subgenotype differences between viral strains that differ in the severity of disease they cause, and possible virulence sites or phenotype-related sites in the genome of CVA6, CVA10 and CVA16.

Materials and methods

Clinical diagnosis, sample collection, and testing by real-time RT-PCR

According to the diagnostic criteria newly defined by the National Health Commission of the People’s of Republic of China in 2018, HFMD manifests with fever and, at the onset, at least one of the following features: maculopapular or vesicular rash on the palms and/or soles and vesicles or ulcers in the mouth. If the disease involves central nervous system symptoms, such as lethargy, ease of being startled, delirium, headache, vomiting, limb tremors, myoclonus, abnormal eye movements, and acute flaccid paralysis, it is diagnosed as severe HFMD. If the serious disease is complicated by at least one of the following signs, it is considered critical HFMD: cold sweats, cold extremities, mottled skin, increased blood pressure, tachycardia, tachypnea, cyanosis, cough with pink foamy or bloody sputum, hypotension, cardiovascular collapse, convulsions, and coma [22].

From 2014 to 2018, a total of 137 fecal specimens from 137 patients with severe or critical HFMD were collected by pediatricians from Shenzhen Children’s Hospital and sent to the Institute of Pathogen Biology at the Shenzhen CDC. Each specimen contained a standardized report that recorded information on patient demographics and clinical findings. Two hundred μl of supernatant from each treated specimen was used to extract viral RNA using a High Pure Viral RNA Kit (Roche) according to the manufacturer’s instructions. Then, the viral RNA was simultaneously screened for EV-A71, CVA16, and pan-enterovirus (PE) by a real-time reverse transcription PCR method (RT-PCR). A known negative specimen was included in each extraction run, and the extract was tested by real-time RT-PCR along with the clinical specimens to monitor for cross-contamination.

Determination of virus based on RT-snPCR and sequencing

For specimens that were negative for both EV-A71 and CVA16 but positive for PE, the strains were typed using two rounds of seminested PCR (snPCR). First, cDNA was synthesized in a 10-μl reaction mixture using a PrimeScript II 1st Strand cDNA Synthesis Kit (Takara, Dalian, China). Second, the primers HEVAS1495 (5’-CAGTCTTCATCAGTTACCCTGGTIATICCITGGATIAGYAAC

AC-3’), HEVAR2807 (5’-CAGTTACACCGGGCAATCGTGTCRCAICCYTGIGCI

GTRG-3’), and HEVAR2C (5’-CGGTGYTTGCTCTTGAACTGCATG-3’), which target the complete VP1 gene, were used to detect EV-A sequences in 5 μl of cDNA solution [23]. Third, two pairs of primers (224 [5’-GCIATGYTIGGIACICAYRT-3’]/222 [5’-CICCIGGIGGIAYRWACAT-3’] and AN89 [5’-CCAGCACTGACAGC

AGYNGARAYNGG-3’]/AN88 [5’-TACTGGACCACCTGGNGGNAYRWACAT-3’]) targeting part of the VP1 gene were used to detect other EV sequences in 5 μl of cDNA solution [24]. Amplified DNA was purified and sequenced by Sangon Biotech (Shanghai, China) Co., Ltd. using the Sanger method.

Virus isolation

Enterovirus-positive specimens were selected for virus isolation. Viral strains were isolated in rhabdomyosarcoma (RD) cell cultures according to the protocol established by the China CDC. Briefly, a mixture of specimen, phosphate-buffered saline (PBS, pH 7.4), benzylpenicillin (100 U/mL), streptomycin sulfate (100 mg/mL), and chloroform (10% [v/v]) was prepared by shaking with glass beads. Two hundred μl of supernatant was prepared after centrifugation at 1,500 × g, and added to 24-well plates with a monolayer of RD cells. The cultures were incubated at 37°C in a 5% CO2 atmosphere. The viral cultures were considered positive if a cytopathic effect (CPE) was observed after three passages.

Complete VP1 gene and near-full-length genome sequencing

The 24 EV-A71 and six CVA16 strains isolated from severe HFMD patients, collected between 2014 and 2017, were selected randomly for amplification of complete VP1 sequences or nearly complete genome sequences. Two hundred μl of supernatant from cell cultures was used to extract viral RNA using a High Pure Viral RNA kit. The complete VP1 sequences of EV-A71 strains were amplified by one-step RT-PCR with the primers EV71-VP1-S (5’-GCAGCCCAAAAGAACTTCAC-3’) and EV71-VP1-A (5’-AAGTCGCGAGAGCTGTCTTC-3’) [25]. The complete VP1 sequences of CVA16 strains were amplified by one-step RT-PCR with the primers CA16-ZF (5’-ATTGGTGCTCCCACTACAGC-3’) and CA16-ZR (5’-GCAAGGTGCCGATTCACTACCCT-3’) [26]. Amplified complete VP1 products were purified and sequenced by Sangon Biotech (Shanghai, China) Co., Ltd. using the Sanger method.

The nearly complete genome sequences of EV-A71, CVA6, CVA16, CVA10 and CVA2 strains were amplified using a pair of universal primers, EVA-F30 (5’-TTAAAACAGCCTGTGGGTTGTACCCACCCA-3’) and EVA-R36 (5’-GCTATTCTGGTTATAACAAATTTA CCCCCACCAGTC-3’), by a one-step RT-PCR method as described previously [27]. Amplified nearly complete genome products were sequenced using a primer-walking method (Takara, Dalian, China). The nearly complete genome sequences were assembled from the generated contigs using the program Sequencher version 4.9.

Phylogeny-based genotyping and genotype comparison between strains associated with different disease severity

Genotyping of EV-A17, CVA6, CVA16, CVA10, and CVA2 was carried out following the methods recommended by China CDC [28,29,30,31,32]. Reference VP1 sequences were obtained from the GenBank database and covered all known genogroups and subgenotypes. Some reference sequences with known disease severity were selected for comparative analysis of subgenotypes of viral strains linked to different disease severity. Phylogenetic trees were constructed using the neighbor-joining method (NJ) implemented in MEGA 7.0.26 [33]. The nucleotide substitution model used for phylogenetic analysis was the p-distance model. Bootstrap analysis (1,000 replicates) was used to evaluate the statistical support for the trees.

Comparative genomic analysis of strains associated with severe and mild HFMD

The information for CVA6, CVA16 and CVA10 strains with known disease background was obtained from previous reports [4, 5, 7, 27, 30, 34,35,36,37,38,39,40] and from the GenBank database. The following nearly complete genome sequences were obtained from GenBank: five CVA6 strains associated with severe HFMD and 28 CVA6 strains associated with mild HFMD; 23 CVA16 strains associated with severe HFMD and 20 CVA16 strains associated with mild HFMD; five CVA10 strains associated with severe HFMD and 26 CVA10 strains associated with mild HFMD. Details about these strains are listed in Supplementary Material Table S1. The 5’ untranslated regions (5’UTRs), deduced amino acid sequences, and 3’ untranslated regions (3’UTRs) of strains associated with severe and mild HFMD were compared. The sequences were analyzed using BioEdit (version 7.2.5) and ClustalX (version 2.0.12).

Statistical analysis

Statistical analysis was performed using Excel 2007 and SPSS 22.0. Statistical differences between proportions were analyzed by the chi-squared test. A P-value less than 0.05 was regarded as statistically significant.

Data availability

Thirty-seven VP1 sequences of eight enterovirus types determined in this study have been deposited in the DDBJ/ENA/GenBank database under the accession numbers MH716145-MH716181. The nearly complete genome sequences of 17 EV strains associated with severe HFMD were deposited in the DDBJ/ENA/GenBank database under the following accession numbers: EV-A71 (KT428644-KT428650), CVA16 (KM215267 and KX595291-KX595294), CVA6 (MH716144 and MN032612), CVA10 (KX595288 and KX595290), and CVA2 (KX595284).

Results

Epidemiology of severe HFMD in Shenzhen, China, from 2014 to 2018

From 2014 to 2018, a total of 137 severe HFMD cases, including 10 fatal cases, were reported to Shenzhen CDC (84 in 2014, 34 in 2015, nine in 2016, eight in 2017, and two in 2018). The number of severe HFMD patients decreased significantly over time. The male/female ratio of severe HFMD patients was 1.9:1(90:47) (Fig. 1A). The majority of severe HFMD patients were aged 0-3 years (74.5%, 102/137). In most of the cases, the patient was 1-2 years old (37.2%, 51/137) (Fig. 1B). Severe HFMD infections (72.3%, 99/137) mainly occurred between April and June, with a peak in May during this 5-year period (Fig. 1C).

Fig. 1
figure 1

Distribution of 137 severe HFMD cases by gender (A), age (B), and month (C) in Shenzhen, China, during 2014-2018

Full size image

Etiology of severe HFMD in Shenzhen, China, from 2014 to 2018

The EV detection rate in severe HFMD cases was 97.6%, 94.1%, 88.9%, 87.5%, and 100% in 2014, 2015, 2016, 2017, and 2018, respectively (Table 1). Overall, the EV detection rate for severe HFMD was 95.6% (131/137) in this 5-year period. A total of eight EV types were identified to be associated with severe HFMD infections during 2014-2018. EV-A71 was the most frequently detected type, with a detection rate of 70.8% (97/137). After EV-A71, the other seven types were CVA6 (8.8%, 12/137), CVA16 (7.3%, 10/137), CVA2 (2.9%, 4/137), CVA10 (2.9%, 4/137), CVA5 (0.7%, 1/137), E6 (0.7%, 1/137), and E18 (0.7%, 1/137). One strain was untyped (0.7%, 1/137). Of the 10 fatal cases, nine patients (90%, 9/10) were positive for EV-A71, and the other one (10%, 1/10) was positive for CVA16 (Table 2).

Table 1 Detection rates and distribution of the eight enterovirus types among 137 patients with severe hand, foot and mouth disease in Shenzhen, China, 2014-2018
Full size table
Table 2 Demographic features and serious complications in 137 patients with severe HFMD associated with the different enterovirus types
Full size table

Demographic features and complications of severe HFMD in Shenzhen, China, 2014-2018

The mean age of EV-A71-infected patients was higher than the mean age of those infected with CVA6 or CVA16 (Table 2). However, the male/female ratio of EV-A71-infected patients was lower than that of patients infected with CVA6 or CVA16. Complications including myoclonic jerk, aseptic encephalitis, tachycardia, tachypnea, lethargy, vomiting, abnormal eye movements, aseptic meningitis, ease of being startled, acute flaccid paralysis (AFP), limb tremors, neurogenic pulmonary edema, myocarditis, and brainstem encephalitis were observed in 137 patients with severe HFMD (Table 2). The most frequent complication in severe HFMD cases with proven EV infection (n = 130) was myoclonic jerk (60%, 78/130). Other common serious complications included aseptic encephalitis (37.7%, 49/130), tachycardia (33.1%, 43/130), and tachypnea (30.8%, 40/130). Abnormal eye movements, ease of being startled, AFP, and neurogenic pulmonary edema were observed only in EV-A71-infected patients.

Relationship between the subgenotype and disease severity

Phylogenetic analysis based on the complete VP1 gene indicated that all of the EV-A71 strains from between 2014 and 2017 belonged to subgenotype C4a (Fig. 2A), which included viral strains from healthy children, mild HFMD patients, severe HFMD patients, and lethal cases, suggesting that there was no association between the subgenotype of EV-A71 and disease severity. All the CVA6 strains from 2014, 2015 and 2018 were classified into subgenotype D3a based on phylogenetic analysis of the VP1gene (Fig. 2B). No specificity was observed between the subgenotypes of CVA6 strains from healthy children, mild HFMD patients, and severe HFMD patients. Phylogenetic analysis based on the complete VP1 gene showed that five CVA16 strains collected in 2014 and 2016 belonged to subgenotype B1b, and another strain collected in 2014 belonged to subgenotype B1a (Fig. 3A). There was no specific association of any of the subgenotypes of CVA16 with mild HFMD, severe HFMD, or lethal infection. Phylogenetic analysis based on the complete VP1 gene indicated that four CVA10 strains from 2014 and 2015 clustered in subgenotype C2 (Fig. 3B), and there was no association between the subgenotype of CVA10 and disease severity. Four CVA2 strains from 2015 and 2016 were assigned to subgenotype D2 based on phylogenetic analysis of the VP1 gene (Fig. 3C), and no specific association was observed between any of the subgenotypes of CVA2 and mild HFMD or severe HFMD.

Fig. 2
figure 2

Neighbor-joining phylogenetic trees for EV-A71 (A) and CVA6 (B) based on complete VP1 sequences. The nucleotide substitution model used was the p-distance model. One thousand bootstrap replicates were used for construction of the phylogenetic trees; values >70% are shown. The scale bar represents a genetic distance of 0.05 nucleotide substitutions per site. The symbol “□” indicates strains from healthy children, “○” indicates strains from patients with mild HFMD, “●” indicates strains from patients with severe HFMD from this study, and “■” indicates strains from fatal cases from this study. The reference sequences are labeled with GenBank accession no./country/year

Full size image
Fig. 3
figure 3

Neighbor-joining phylogenetic trees based on complete VP1 sequences for CVA16 (A), CVA10 (B) and CVA2 (C). The p-distance model was used to construct the trees. Bootstrap analysis (1,000 replicates) was used for statistical support of the trees, and only bootstrap support values > 70% are shown. The scale bar represents a genetic distance of 0.05 nucleotide substitutions per site. The symbol “□” indicates strains from healthy children, “○” indicates strains from patients with mild HFMD, “●” indicates strains from patients with severe HFMD from this study and from the GenBank database, and “■” indicates lethal strains from this study. The reference sequences are labeled with GenBank accession no./country/year

Full size image

Origin of EV strains from this study

In a BLAST search, the nearly complete genome sequences of seven EV-A71 strains from this study showed the closest genetic relationship to those of EV-A71 strains from mainland China sharing 97.02%-99.89% nucleotide (nt) sequence identity with these strains (Table S2). The other 17 complete VP1 sequences of EV-A71 strains showed 98.32% to 99.78% nt sequence identity to those of the most closely related strains from mainland China. The nearly full-length genome sequence of the strain CVA6/sHFMD14/Shenzhen/2015 showed 98.04% nt sequence identity to that of the closest strain, CVA6/S2727/BJ/CHN/2014, from mainland China, and the strain CVA6/sHFMD02/Shenzhen/2018 showed 97.58% nt sequence identity to the closest strain, SHAPHC6069/SH/CHN/15, from mainland China (Table S2). The other 10 complete VP1 sequences of CVA6 strains and their genetically closest sequences from mainland China shared 99.34% to 100% nt sequence identity. The nearly complete genome sequences of five CVA16 strains from this study showed the closest genetic relationship to those of CVA16 strains from mainland China and shared 97.65%-98.15% identity with the closest strains (Table S2). The complete VP1 sequence of the strain CVA16/sHFMD03/Shenzhen/2016 showed the closest genetic relationship to that of the strain CVA16/Shenzhen87/CHN/2016 from Shenzhen, with 100% identity. The nearly complete genome sequence of the strain CVA2/Shenzhen21/CHN/2015 showed the closest genetic relationship to that of the strain BJ13-53/BJ/CHN/2013 from Beijing, with 97.76% identity. The other three complete VP1 sequences of CVA2 strains showed the closest genetic relationship to those of CVA2 strains from mainland China and Hong Kong, with 97.63%-98.53% identity. The nearly complete genome sequences of the strains CVA10/Shenzhen18/CHN/2014 and CVA10/Shenzhen10/CHN/2015 shared 99.16% and 98.11% identity with those of the closest strains CV-A10/P911/2013/China and CV-A10/P670/2013/China, respectively, from mainland China (Table S2). The other two complete VP1 sequences of CVA10 strains showed 99.66% identity to those of the closest strains from mainland China. The one CVA5 strain from 2015 showed 98.03% nt sequence identity to both CVA5 strains associated with an epidemic outbreak in Taizhou, China, in 2015 and the strain JB141330144-CA5 collected in Shenzhen, in 2013, suggesting that these CVA5 strains collected in 2015 might have originated from the 2013 Shenzhen strain JB141330144-CA5 (Table S2). The partial VP1 sequence of an E6 strain from the study shows the closest genetic relationship to that of the 2012 Beijing strain FT2012HFM-155 from environmental sewage, with 94.53% identity (Table S2). The partial VP1 sequence of an E18 strain from this study showed the closest genetic relationship to that of the Hunan strain S1871B from a patient with HFMD, with 98.19% identity (Table S2).

Analysis of virulence- or phenotype-associated sites in CVA6, CVA16, and CVA10

The deduced amino acid (aa) sequences of the polyproteins of the five CVA6 strains associated with severe HFMD and the 28 CVA6 strains associated with mild HFMD were aligned and analyzed. Sites where amino acids differ in more than two strains are defined as variant sites [41]. The proportions of variant residues in six such sites display significant differences between five CVA6 strains associated with severe HFMD and 28 CVA6 strains associated with mild HFMD (Table 3). At position 2 of the 2A protease (2Apro), the proportion of R/K in CVA6 strains associated with severe HFMD is different from that in CVA6 strains associated with mild HFMD (χ2=4.1026, p=0.0428). The residue R is more likely to appear in CVA6 strains associated with severe HFMD than in CVA6 strains associated with mild HFMD. The proportions of R/K and N/S are higher in severe CVA6 strains than in mild CVA6 strains at the positions 41 and 48, respectively, of the 2C protein (χ2 = 7.9974, p = 0.0047 and χ2 = 7.0714, p = 0.0078, respectively). The constituent ratios of residues A and V are higher in CVA6 strains associated with severe HFMD than in CVA6 strains associated with mild HFMD at positions 15 and 68, respectively, of the 3A protein (χ2 = 8.2598, p = 0.0041 and χ2 = 4.1026, p = 0.0428, respectively). At position 67 of the 3D polymerase (3Dpol), the proportion of I/L is higher in CVA6 strains associated with severe HFMD than in CVA6 strains associated with mild HFMD (χ2 = 5.3049, p = 0.0213). In the 5’ untranslated region (5’UTR), the 33 sequences had same length at nucleotide positions 183-747. Among these 565 nucleotide positions, bases in four positions display significant differences between CVA6 strains associated with severe HFMD and CVA6 strains associated with mild HFMD. The proportions of U/C, C/U, A/C and U/C are higher in CVA6 strains associated with severe HFMD than in CVA6 strains associated with mild HFMD at positions 306, 502, 650 and 707, respectively, in the 5’UTR (χ2 = 12.6818, p = 0.0004; χ2 = 12.6818, p = 0.0004; χ2 = 9.9754, p = 0.0016 and χ2 = 9.9754, p = 0.0016, respectively) (Table 3). In the 3’UTR, the 33 sequences had same length at nucleotide positions 1-38. Among the 38 nucleotide positions of the 3’UTR, bases at two positions display significant differences between CVA6 strains associated with severe HFMD and CVA6 strains associated with mild HFMD. The proportions of C/U and C/U are higher in CVA6 strains associated with severe HFMD than in CVA6 strains associated with mild HFMD at positions 7 and 27, respectively, in the 3’UTR (χ2 = 8.2598, p = 0.0041 and χ2 = 9.0659, p = 0.0026, respectively) (Table 3). The above variant sites are potentially associated with the virulence or clinical phenotype of CVA6. When sites 650 and 707 in the 5’UTR, sites 41 and 48 in the 2C protein, site 15 in the 3A protein, and sites 7 and 27 in the 3’UTR were used as virulence or phenotype sites to identify CVA6 strains associated with severe HFMD, four (80%) of the five CVA6 strains associated with severe HFMD were uniformly identified (Table 3).

Table 3 Possible virulence- or phenotype-associated sites revealed by comparative genomic analysis of five CVA6 strains associated with severe HFMD and 28 CVA6 strains associated with mild HFMD
Full size table

Comparative genomic analysis of the 23 CVA16 strains associated with severe HFMD and the 20 CVA16 strains associated with mild HFMD revealed three significant differences in the 5’UTR, four in the entire coding region, and one in the 3’UTR (Table 4). In the 5’UTR, the proportions of C/U, G/A and G/A are higher in CVA16 strains associated with severe HFMD than that in CVA16 strains associated with mild HFMD at positions 200, 227 and 747, respectively (χ2 = 4.7685, p = 0.0290; χ2 = 4.7685, p = 0.0290; and χ2 = 5.9595, p = 0.0146, respectively). At position 91 of the 2B protein, the proportion of F/L is higher in CVA16 strains associated with severe HFMD than in CVA16 strains associated with mild HFMD (χ2 = 6.0635, p = 0.0138). At position 11 of the 3A protein, the proportion of N/S is higher in CVA16 strains associated with severe HFMD than in CVA16 strains associated with mild HFMD (χ2 = 5.7336, p = 0.0166). The proportions of A/V and I/T are higher in CVA16 strains associated with severe HFMD than in CVA16 strains associated with mild HFMD at positions 74 and 92, respectively, of 3Dpol (χ2 = 4.9199, p = 0.0265 and χ2 = 8.3198, p = 0.0039, respectively). The constituent ratio of the base U is higher in CVA16 strains associated with severe HFMD than in CVA16 strains associated with mild HFMD at position 30 in the 3’UTR (χ2 = 4.9199, p = 0.0265).

Table 4 Possible virulence- or phenotype-associated sites revealed by comparative genomic analysis of 23 CVA16 strains associated with severe HFMD and 20 CVA16 strains associated with mild HFMD
Full size table

A total of 10 significant variant sites were found in the UTR and polyprotein region by comparative genomic analysis of the five CVA10 strains associated with severe HFMD and the 26 CVA10 strains associated with mild HFMD (Table 5). The proportions of A/G, U/C, G/A and U/C are higher in CVA10 strains associated with severe HFMD than in CVA10 strains associated with mild HFMD at positions 220, 423, 691 and 731, respectively, of 5’UTR (χ2 = 7.5164, p = 0.0061). The proportions of V/A and V/I in CVA10 strains associated with severe HFMD are higher than in CVA10 strains associated with mild HFMD at positions 23 and 283, respectively, of VP1 (χ2 = 7.5164, p = 0.0061 and χ2 = 7.5164, p = 0.0061, respectively). The proportion of V/I is higher in CVA10 strains associated with severe HFMD than in CVA10 strains associated with mild HFMD at position 8 of the 3A protein (χ2 = 7.5164, p = 0.0061). The proportions of V/I, G/D and H/R in CVA10 strains associated with severe HFMD are higher than in CVA10 strains associated with mild HFMD at positions 15, 32 and 33, respectively of the 3C protease (3Cpro) (χ2 = 7.5164, p = 0.0061; χ2 = 11.2432, p = 0.0008; and χ2 = 7.5164, p = 0.0061, respectively). No significant virulence-associated site was found in the 3’UTR of CVA10. When all of the above variant sites were used as virulence or phenotype markers to identify CVA10 strains associated with severe HFMD, four (80%) of the five CVA10 strains associated with severe HFMD were uniformly identified.

Table 5 Possible virulence- or phenotype-associated sites revealed by comparative genomic analysis of five CVA10 strains associated with severe HFMD and 26 CVA10 strains associated with mild HFMD
Full size table

Discussion

There were 93, 110, 185, 52, 15, 84, 34, 9, 8 and 2 cases of severe HFMD reported to Shenzhen CDC in 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017 and 2018, respectively [3]. Peaks of severe HFMD infections occurred in 2011 and 2014 in Shenzhen. The number of cases of severe HFMD decreased significantly after 2016 in Shenzhen. There are several possible explanations for this. First, an EV-A71 vaccine was introduced in mainland China in 2016. Second, emerging CVA6 and reemerging CVA16 strains have become the main pathogens associated with mild HFMD in recent years in Shenzhen, while the prevalence of EV-A71 is very low (data not shown). Third, public awareness of HFMD has increased. Thus the disease can be discovered and treated at an early stage. This study revealed that the number of male patients is higher than that of female patients, the majority of patients with severe disease were under the age of three, and the peak of severe HFMD infections occurred in May. These epidemiological features were consistent with those reported for patients with mild HFMD in southern China [1, 30].

A previous study revealed that CVA2, CVA4, CVA5, CVA6, CVA10, CVA16, EV-A71, CVB2, CVB4, CVB5 and E9 were associated with severe HFMD in Shenzhen from 2009 to 2013 [3]. The predominant pathogen was EV-A71, followed by CVA6 and CVA16 during this period. The current study indicates that the distribution of the main causative agents of severe HFMD in 2014-2018 was similar to that in 2009-2013. No serious cases related to CVB infection were found during 2014-2018 in Shenzhen, but two new EV types, E6 and E18, were found to be associated with serious HFMD infections. Of 592 patients with severe HFMD reported between 2009 and 2018 [3], 72.3% were infected with EV-A71, 6.3% exhibited CVA6 infections, and 4.9% had CVA16 infections (Table 6). The detection rate of EV-A71 was significantly higher than that of CVA6 and CVA16 in patients with severe HFMD (χ2 = 869.6442, p < 0.0001). Although EV-A71 was the main pathogen associated with severe HFMD, all of the severe infections caused by it were sporadic cases. The severe HFMD infections caused by other EV types were also sporadic cases, but the number of severe HFMD cases associated with CVA6 and CVA10 infections has increased in recent years in mainland China [5, 8,9,10]. Therefore, it is necessary to strengthen surveillance of severe HFMD to prevent outbreaks.

Table 6 Comparison of detection rates and serious complications of EV-A71, CVA6 and CVA16 in 592 patients with severe hand, foot and mouth disease in Shenzhen, China, 2009-2018
Full size table

More than 10 complications were observed in 494 severe HFMD patients with proven EV infections (Table 6). The most frequent complication was myoclonic jerk (61.9%, 306/494), and other common severe complications included aseptic encephalitis (47.0%, 232/494), tachypnea (29.8%, 147/494), vomiting (28.1%, 139/494), tachycardia (24.5%, 121/494), and lethargy (17.0%, 84/494). Gonzalez et al. reported that myoclonic jerk is a significant indicator of central nervous system complications in EV-A71-infected HFMD patients [17], and this was also found in our study. Furthermore, our study showed that myoclonic jerk is also the most frequent complication in severe CVA6 and CVA16-associated HFMD. There was no significant difference between the frequencies of serious complications in EV-A71-infected, CVA6-infected, and CVA16-infected patients, except for those of vomiting and abnormal eye movements, which were more common in severe EV-A71-associated disease than in severe CVA6- or CVA16-associated disease, suggesting that vomiting and abnormal eye movements could be used as clinical parameters to distinguish severe EV-A71 and CVA6 or CVA16 serious infections (Table 6).

Sequence analysis demonstrated that all of the strains included in this study showed the closest genetic relationship to strains from mainland China. Hence, the EV strains associated with severe HFMD circulating in Shenzhen are considered indigenous strains. Phylogenetic analysis indicated that there is no relationship between the subgenotype of the virus and disease severity. A previous study revealed no association between disease severity and patient gender or age [6]. However, EV displays considerable phenotypic variation. Genome sequence comparisons showed that four variable sites in the 5’UTR, six in the polyprotein region, and two in the 3’UTR are likely to be associated with the virulence or clinical phenotype of CVA6. The four variable sites are located in two functional regions of the 5’UTR: the type I internal ribosome entry site (IRES) element (U306C and C502U) and the hypervariable region (A650C and U707C), which regulate the initiation of viral protein translation [41, 42]. Several studies have shown that variations in multiple sites in the IRES (C158U, 272 and 488) and the hypervariable region (700) of EV-A71 influence its virulence [19, 20]. One variable site (R2K) is located in 2Apro, which performs the posttranslational proteolytic processing of the viral polyprotein and also cleaves several host-cell proteins, which promotes the production of new virus particles and allows the virus to evade cellular antiviral immune responses [43,44,45]. The residue K at position 68 of 2Apro is more frequently found in severe cases of EV-A71 [41]. Two variable sites (R41K and N48S) are located in the 2C protein, which has ATPase-dependent helicase and ATPase-independent RNA chaperone activity [44, 45], but no virulence-associated sites have been reported in 2C of EV-A71. Two variable sites (A15T/M/K and V68I) are located in 3A, which interacts with acyl-CoA-binding domain-containing protein 3 (ACBD3) to recruit the lipid kinase phosphatidylinositol 4-kinase-β (PI4KB) to genome replication sites [44, 45]. There is no virulence-associated site reported in 3A of EV-A71. One variable site (I67L) is located in 3Dpol, which functions as a RNA-dependent RNA polymerase (RdRp) and plays an important role in changing viral virulence [19, 44, 45]. The I251T variation in the 3Dpol of EV-A71 affects viral virulence and temperature sensitivity [20]. Two variable sites (C7U and C27U) are located in 3’UTR, which is involved in negative-strand RNA synthesis and viral translation [46]. Nevertheless, no virulence-associated sites were found in the 3’UTR of EV-A71. When seven variable sites (sites 650 and 707 in the 5’UTR, sites 41 and 48 in 2C, site 15 in 3A, and sites 7 and 27 of 3’UTR) were used in combination as virulence markers, it was possible to identify 80% of CVA6 strains associated with severe HFMD. Four CVA6 strains associated with severe HFMD exhibited consistent differences in seven of 10 variable sites compared with another CVA6 strain associated with severe HFMD. These seven sites probably act synergistically to influence the virulence or clinical phenotype of CVA6. Comparative analysis revealed eight significant sequence differences between 23 CVA16 associated with severe HFMD and 20 CVA16 strains associated with mild HFMD. The two variable sites (C200U and G227A) in the 5’UTR were present in the same nine strains associated with severe HFMD. These two sites might act synergistically to affect the virulence or clinical phenotype of CVA16. This pattern was not observed at other variable sites. One variable site is located in 2B, which is an ion channel protein with two transmembrane domains (TM1 and TM2) [44, 45]. The 2B protein has viroporin activity and contributes to lipid kinase phosphatidylinositol 4- kinase- β (PI4KB) recruitment, but its exact function is still unknown. No virulence-associated sites have been reported in 2B of EV-A71. Ten significant differences were observed between five CVA10 associated with severe HFMD and 26 CVA10 strains associated with mild HFMD. Interestingly, four severe CVA10 strains exhibit consistent difference in the above 10 sites compared to the other CVA10 strain associated with severe HFMD. This suggests that these 10 sites probably have a synergistic effect in influencing the virulence or clinical phenotype of CVA10. Four variable sites are located in two functional regions of the 5’UTR: the IRES element (A220G and U423C) and the hypervariable region (G691A and U731C). Two variable sites (V23A and V283I) are located in capsid protein VP1, which determines the antigenicity of the virus and binds to cell receptors and mediates entry into cells [44, 45]. Previous studies have shown that multiple sites or residues in VP1 of EV-A71, including 31G, 31D, L79R, 98K/E, 104D, 107A, Q145E/G/R, 164E, and K244E, are associated with its virulence [20]. Three sites (V15I, G32/D and H33R) are located in 3Cpro, which is a cysteine protease and functions like 2Apro [43,44,45]. Studies have shown that the variation N69D in 3Cpro is associated with the virulence of EV-A71, and variations at residue 79 (T79V/I/A) of 3Cpro are associated with the clinical severity and replication efficiency of EV-A71 [47, 48]. Individual variations can therefore potentially change the virulence or clinical phenotype of EV-A71 by affecting viral replication, translation, temperature sensitivity, and interaction with host cells. Hence, we presume that variations in genomes of CVA6, CVA16 and CVA10 influence their virulence or clinical phenotype in the same way they affect EV-A71.

The main limitation of the comparative genomic analysis of this study is the small number of viral genome sequences compared, especially those of strains associated with severe HFMD. In spite of this, statistical analysis of differences between strains associated with severe HFMD and strains associated with mild HFMD was still performed in the present study. Virulence site analysis for other EV types was not conducted due to the limited number of available genome sequences with known disease background. In this study, we only explored the possible relationship between the virus genome and disease severity. Reverse genetic experiments are needed to further identify potential virulence sites. The clinical phenotype resulting from viral infections is determined by complex interactions of multiple factors [49]. Besides the virus, host factors, the intestinal virome, coinfection with more than two EV types, and the environment can also influence HFMD outcome [6, 50,51,52,53,54,55]. No coinfection was found in patients with severe HFMD based on the traditional experimental methods used in this study. Further studies are needed to help us understand the mechanism of severe HFMD.

In summary, EV-A71 was the main pathogen associated with severe HFMD in Shenzhen, China, during 2014-2017, but various other EV types can also cause serious complications. The EV strains associated with severe HFMD circulating in Shenzhen in past years appear to be indigenous strains. No association between the subgenotype of the virus and disease severity was found. However, virulence or phenotype-associated sites potentially exist in the genomes of CVA6, CVA16 and CVA10.