Introduction

Coxsackievirus A16 (CVA16) is a small, non-enveloped virus containing a single-stranded positive-sense RNA genome. It is one of the major pathogens associated with hand, foot, and mouth disease (HFMD) in infants and young children [1,2,3]. First isolated in South Africa in 1951 [4], CVA16 infection-caused outbreaks of HFMD have been reported worldwide and have been particularly widespread in the Asia-Pacific region [5,6,7,8,9]. The majority of HFMD with CVA16 infections present with mild symptoms [5, 10]. CVA16 infection may contribute to severe complications, including cardiac and pericardial disease, pulmonary complications, spontaneous abortion, myocarditis and pneumonia, aseptic meningitis, rhombencephalitis, and death [1, 5, 11,12,13,14]. Previous studies found CVA16 infection to be responsible for approximately 21% of severe HFMD cases with neurological complications [15]. Viral infection by CVA16 is a serious public health issue in affected areas.

The first HFMD case was reported in Shanghai in 1981. Since then, HFMD cases have been reported in many provinces and municipalities in China [1, 6, 9, 16,17,18,19]. In mainland China, approximately 50% of all laboratory-confirmed enterovirus infections in HFMD cases are caused by CVA16 [6]. In 2008, HFMD was classified as a class C infectious disease in mainland China [3, 20]. The frequent recombination and mutation of virus genes has led to HFMD outbreaks associated with CVA16 presenting with different symptoms in different regions [6, 21]. There is currently no CVA16 antiviral treatment available. Vaccination is considered the most effective method to control the transmission and infection rate of viruses. However, there is no CVA16 vaccine available for public use [22], and the listed EV71 vaccine has no significant protective effect for HFMD caused by CVA16 [23]. The CVA16 genome contains approximately 7410 nucleotides, including a single open-reading-frame (ORF) flanked by 5/3-untranslated regions (UTR). The ORF has 6579 nucleotides and encodes a polyprotein of 2193 amino acids, which is composed of three protein precursors: P1, P2, and P3. The P1 polyprotein precursor is processed into four structural proteins (VP1, VP2, VP3, and VP4), while P2 and P3 are precursors of the seven nonstructural proteins (P2–2A, 2B, 2C; P3–3A, 3B, 3C, and 3D) [24, 25]. VP1, VP2, and VP3 capsid proteins are exposed on the viral surface, while VP4 is usually situated on the inner portion. There are many B-lymphocyte antigenic epitopes and virulence-associated sites in the viral structural proteins [26,27,28,29,30]. Investigations into the geographical distribution and genetic characteristics of the gene in the P1 coding region of CVA16 associated with HFMD in China may provide basic information for the prevention and control of HFMD.

Materials and methods

Objectives

To assess the geographical distribution and genetic characteristics of the gene in the P1 coding region of CVA16 associated with HFMD in China, we performed an extensive genetic analysis using all available P1 coding region sequences of CVA16 in the public database (dated December 2017 or earlier).

Study design

Sequence collection and preliminary treatment

All CVA16 P1 coding region sequences associated with HFMD in China with known collection dates and isolate regions, and 36 complete VP1 gene sequences of CVA16 strains used as references, were obtained from the GenBank public database on the National Center for Biotechnology Information (NCBI) PubMed website (http://www.ncbi.nlm.nih.gov/genbank/). The accession numbers and specific sequence data are included in Supplementary Table 1. Molecular Evolutionary Genetics Analysis software (MEGA; ver. 6.06) was used to edit and align the nucleotide and deduced amino acid sequences with ClustalW. Repeat sequences, defined as identical in terms of sequence similarity, isolation time, and isolation location, were removed from the dataset to generate a non-redundant dataset.

Phylogenetic analysis, substitution rates and recombination analysis

The phylogeny of the VP1 genes of CVA16 was examined with neighbor-joining (N-J) trees reconstructed with MEGA [6], using Kimura’s 2-parameter substitution model with 1000 bootstrap replicates. To compare the rates between synonymous (Ka) and nonsynonymous (Ks) substitutions, we estimated Ka and Ks using alignments of the P1 coding genes in MEGA and analyzed them in DnaSp [31]. To assess the influence of recombination on the CVA16 diversity, we analyzed the P1 coding region sequences with the program RDP [32].

Analysis of the homology identity of nucleotide and deduced amino acid sequences

BioEdit (ver. 7.0.9.0; http://www.mbio.ncsu.edu/BioEdit/bioedit.html) was used to perform the homology identity analysis of the nucleotides and deduced amino acid sequences of CVA16 strains in different gene regions of the P1 coding region.

Entropy-based analysis of the diversity of the structural protein-encoding regions (P1) and amino acid mutations analyses

The entropy value of each amino acid site on the structural protein-encoding regions (P1) of Chinese CVA16 strains were calculated with BioEdit (ver. 7.1.9.0). MEGA was used to analyze the deduced amino acid sequence variation sites in the P1 region of the CVA16 strains.

Results

Genotyping based on the complete VP1 gene sequences of CVA16

In total, 184 sequences with a complete VP1 region of CVA16 were analyzed, including 148 Chinese CVA16 strains retrieved from the public database dated prior to December 2017 (Supplementary Table S1). The phylogenetic dendrogram illustrates the division of all CVA16 strains into three distinct genogroups: A, B, and D (Fig. 1). Genogroup A included the CVA16 prototype strain and two Chinese strains (EU812514-AH-2008 and KX595295-GD-2014). A comparison indicated a 99.4–100% nucleic acid identity of the complete VP1 gene sequences. The nucleotide difference of the other CVA16 strains relative to the CVA16 prototype strain was greater than 22.8%. The B2 subgenogroup included five Japanese strains and two Malaysian strains isolated between 1981 and 2000. A total of 170 strains belonging to genotype B1 were clustered together and categorized into three subtypes: B1a, B1b, and B1c. Compared with the CVA16 prototype strain, the nucleotide identity of the VP1 gene of all the Chinese CVA16 strains was 75.1–100%. Of the 148 Chinese CVA16 strains isolated from the period 1998–2014, the nucleotide identity of the VP1 gene was 74.8–100%. The majority were grouped in the B1 subgenogroup: 98 strains clustered together with the reference sequences of subgenotype B1b, and 48 strains clustered with the subgenotype B1a reference strains. Two strains, isolated in Anhui in 2008 and Guangdong in 2014, clustered together with the CVA16 prototype strain belonging to genotype A. There were no B1c, B2, or D subgenotypes identified in China (Fig. 1).

Fig. 1
figure 1

Phylogenetic tree constructed using the VP1 nucleotide sequences (891 bp) of Coxsackievirus A16 (CVA16) analyzed by the neighbor-joining (N-J) method with 1000 bootstrap replicates. CVA16 reference sequences are indicated by a red diamond; the EV71 prototype strain (EV71 BrCr) used as the outgroup is shown as a black diamond. Yellow branches indicate genotype A strains; red branches indicate genotype B1b strains; blue branches indicate genotype B1a strains; pink branches indicate genotype B1c strains; green branches indicate genotype B2 strains; and sepia branches indicate genotype D strains. Each strain was named according to the GenBank accession number, country or area of isolation, and year of isolation. USA United States of America, RSA South Africa, MAS Malaysia, JPN Japan, VNM Vietnam, THA Thailand, KOR Korea, AUS Australia, CHN China. TW Taiwan, GD Guangdong, FJ Fujian, AH Anhui, YN Yunnan, BJ Beijing, ZJ Zhejiang, JS Jiangsu, HeN Henan, SD Shandong, SH Shanghai, JL Jilin, Tibet = Xizang; GX Guangxi, HuB Hubei, HeN Henan, HuN Hunan, QH Qinghai, LN Liaoning. (Color figure online)

Full size image

Geographic distribution and genetic recombinant of the CVA16 strains isolated in China

The sampling locations for the CVA16 strains included in the phylogenetic analysis were distributed throughout 18 provinces and municipalities in China (Supplementary Fig. S1). Two strains belonging to the A genotype were identified in Anhui and Guangdong. The majority of the 98 B1b strains and 48 B1a strains were identified in coastal areas. Of the Chinese CVA16 strains, four strains might show inter- and intratypic recombination in the P1 region at different sites: KU163608-HuB-2010, JX068829-HeB-2010, JX986742-HuB-2011, and AY895099-GD-2002 (Table 1 and Supplementary Fig. S2).

Table 1 Recombination analysis of the P1 coding region sequences of Chinese CVA16 stains with RDP4
Full size table

Nucleotide and amino acid identities of the structural protein encoding regions (P1) among the Chinese CVA16 strains

At the P1 gene level, the nucleotide identity of all Chinese CVA16 strains was 73–100%, and the nucleotide divergence between the Chinese CVA16 strains and the CVA16 prototype strain was less than 26.8%. Compared with the CVA16 prototype strain, there were no insertions observed in the P1 encoding regions, although some deletions were observed at the position 2383–2385 and 2580–2586 in the P1 encoding regions of the strain (JX507808-BJ-2012 and AY821796-GD-2004). Among the Chinese CVA16 strains, the nucleotide identities of the VP1, VP2, VP3, and VP4 structural proteins were 74.8–100%, 69.5–100%, 74.3–100%, 65.7–100%, respectively, while the amino acid identities were 88.2–100%, 83.0–100%, 87.6–100%, and 72.4–100%, respectively. In the structural protein-encoding regions, VP1, VP2, VP3, and VP4, nucleotide sequence identities within subgenotype B1b and B1a were greater than 88.3%, 70.7%, 75.3%, and 65.7%, respectively, while the amino acid similarities based on each gene region were greater than 95.2%, 83%, 87.6%, and 72.4%, respectively (Table 2).

Table 2 Comparison of the nucleotide and amino acid identities of the Chinese CVA16 strains (%)
Full size table

Entropy value of each amino acid site on structural protein VP1–VP4

Entropy-based analysis of the diversity of the structural protein-encoding region (P1) sequences with BioEdit indicated that the amino acid sequence of the VP1–VP4 protein encoded by the P1 gene was highly conserved (entropy value < 0.6), although the position aa-295 in the P1 encoding region was a mutable site (entropy = 0.72523) (Fig. 2).

Fig. 2
figure 2

Entropy value of each amino acid site on the VP1–VP4 structural proteins. A higher entropy value or longitudinal axis peak indicates a higher frequency of amino acid mutation sites. A high mutation site is defined as entropy > 1; a mutable site is defined as 0.6 ≤ entropy < 1.0; a conservative site is defined as entropy < 0.6

Full size image

Comparison of substitution rates and amino acid sequence alignment of the structural protein encoding regions (VP1–VP4) of CVA16 strains isolated in China

Nucleotide substitutions were identified throughout the P1 encoding region gene of the CVA16 strains.Comparison of the synonymous (Ks), nonsynonymous (Ka), and Ka/Ks substitution rates between CVA16 BrCr and Chinese CVA16 stains gave values of 0.30176, 0.06350, and 0.2104, respectively. Among the 148 CVA16 strains isolated in China, 88, 76, 58 and 29 amino acid substitutions were identified in the VP1, VP2, VP3, and VP4 sequences, respectively (Table 3). Few amino acid substitutions were observed in the P1 protein-encoding region of CVA16 strains. However, amino acid changes were observed with high frequency in the P1 protein-encoding regions (VP1-T11A (10%), N14S (14%), L23M/V (11%), T98M (16%), V107A (14%), N102D (6.1%), E145V (8.8%), N218D (10%), E241K (22%), T248A/I (6.8%); VP2-I217V (22%), T226A (38%); VP3-N141S/G (5.4%), N240D (15%)) (Supplementary Table S2A–B). Some amino acid changes differ between the Chinese B1a and B1b stains (Table 4).

Table 3 Characteristics of the structural protein encoded by the P1 coding region in Chinese CVA16 strains
Full size table
Table 4 Comparision of the amino acid substitution in the Chinese B1a and B1b stains
Full size table

Discussion

Previous studies have not performed a systematic epidemiologic analysis of CVA16. All complete P1 sequences from the GenBank database were analyzed to identify the genetic characteristics and distribution of CVA16 genotypes in China. Phylogenetic analysis based on VP1 clustered 148 Chinese CVA16 strains into two genotypes: A and B. Two strains belonging to genotype A were isolated in Anhui (2008) and Guangdong (2014). There were 98 strains classified as subgenotype B1b and 48 classified as subgenotype B1a. The B2 and B1c subgenotypes were not identified in China, which is consistent with the findings of previous studies [4, 19, 33]. The B1a and B1b subgenotypes were identified in most Chinese provinces, particularly in the developed coastal regions. The CVA16 genotype strains isolated in China were similar to those identified in neighboring countries, such as Japan, Vietnam, and Thailand, indicating that the Chinese CVA16 strains associated with HFMD outbreaks resulted from a coincident circulation with neighboring countries. The B1c, B2, and D subgenotypes were not detected in China, indicating there may be an important transmission risk for these adventive subgenotypes. Continuing surveillance for new CVA16 genotypes and subgenotypes should be maintained. Previous analyses demonstrated that the CA16 strains are progeny of the CA16 prototype strain and other enterovirus A strains, which showed inter- and intra-typic recombination in the 5′-untranslated region and P2 and P3 non-structural protein encoding regions, while it was rare in the P1 region. This probably occurred because the structural proteins are important for virus replication and correct viral structural assembly [6, 9]. In our study, we determined that some strains (KU163608-HuB-2010, JX068829-HeB-2010, JX986742-HuB-2011, and AY895099-GD-2002) might show inter- and intra-typic recombination in the P1 region. While recombination in the P1 gene region is most often deleterious and likely fatal to the virus [34], the virulence and pathogenicity of these strains might have changed considerably, which may have led to the outbreak of HFMD in China.

Entropy-based analysis of the diversity of the P1 protein-encoding region sequences with BioEdit indicated that the amino acid sequence of the VP1–VP4 protein encoded by the P1 gene was highly conserved among Chinese CVA16 strains.

The mean evolutionary rate of CVA16 is slightly lower than EV71 [4, 35]. As an RNA virus, CVA16 generally has a high mutation rate during viral replication due to the deficiency of proof reading activity [36]. Synonymous (Ks) and nonsynonymous substitution (Ka) patterns are valuable in gene evolution studies. Our results demonstrate that the substitution rate (Ka/Ks) was less than 1, suggesting that most genes are under purifying selection. Among Chinese CVA16 strains, 251 amino acid mutation sites were observed in viral structural proteins encoded by the P1 gene, using the vaccine candidate strain (JN590244-JL-2010) as a reference. These variations may account for the difference in pathogenicity of theses CVA16 strains. The VP1–VP3 proteins are located on the external surface of the capsid and have many functional regions related to immunity, virulence, and pathogenicity [27, 29, 37]. Previous studies have reported many B lymphocyte antigenic epitopes and virulence-associated sites on the viral structural proteins [26, 28, 29, 38]. Amino acid sequence alignment with CVA16 strains indicated that most Chinese CVA16 strains were highly conserved, with few amino acid substitutions in the structural proteins, but several amino acid changes occurred with high frequency (VP1-T11A, N14S, L23M/V, T98M, V107A, N102D, E145V, N218D, E241K, T248A/I; VP2-I217V, T226A; VP3-N141S/G, N240D; VP4-I61M), especially in Chinese B1b stains. These amino acid changes may be associated with CVA16-related HFMD outbreaks in China, necessitating further study of amino acid mutations in CVA16 viral structural proteins. A previous study found six well identified, highly conservative linear neutralizing epitopes (aa 94–108, aa 109–123, aa 163–177, aa 187–201, aa 211–225, aa 271–285) in the VP1 protein of CVA16 [27]. Chinese CVA16 strains exhibited few amino acid changes in these regions, but a high mutation frequency was observed for amino acids at the 98, 107, and 218 positions in the VP1 protein. The conserved peptide segment (TMPTTGTGTGTGNTDGYVN) was at 94–108 in the VP1 protein, which could elicit the production of neutralizing antibodies against CVA16 infection [27]. However, eight CVA16 strains (KC695830-SH-2013, KJ850274-TW-2007, JX127277-FJ-2009, JX127273-FJ-2009, KF193625-GX-2010, JX127259-FJ-2010, JX127260-FJ-2010, and KF193626-FJ-2010) had one amino acid change (N102D) in this neutralizing epitope, mutating the neutral amino acid into an acid amino acid, which may decrease affinity and immunodominance and further affect the combination of neutralizing antibodies or vaccines. Most amino acid changes at a single point were observed on the reported B-lymphocyte epitopes in the VP2 structural protein, except for positions 217 and 226 in VP2; 33 strains had I217V and 55 strains had T226A at those positions, respectively, which may affect their immunogenicity. A previous study identified five B-lymphocyte epitopes (aa4–11, aa56–66, aa87–97, aa139–148, aa176–190, and aa229–237) located in the VP3 structural protein [29, 39]. The aa176–190 of VP3 has been identified as a linear immunodominant epitope that can induce neutralizing antibodies in mice, but not in rabbits [39]. Several amino acid changes were observed in this region, which may lead to decreases in affinity and immunogenicity, further affecting the combination with neutralizing antibodies. Among Chinese CVA16 strains, 29 amino acid site mutations were observed in the VP4 protein, and 14 CVA16 strains had few amino acid mutations in only one reported B lymphocyte epitope (aa 40–54) in the VP4 protein [29]. This indicates that most Chinese CVA16 strains were highly conserved in the VP4, which may be due to the VP4 location within the capsid and less exposure to immune pressure [37].

Multiple sites related to virulence and pathogenicity have been reported on VP1 and VP2. Mutations R51K, K52R, M98T, N102D, T103A, N218D, E241K, T248A/I, V251L and T/H295A in the VP1 protein can decrease the virulence and pathogenicity of CVA16 in mice [23, 40]. Amino acid changes to N102D or E241K in the VP1 decreased virulence, while it changes to M220L increased virulence and the fatality rate in rats [26, 38]. In this study, we observed a high mutation frequency of amino acids at the 98, 107, and 218 positions in the VP1 protein. VP1-L220, located within the neutralization epitope on the GH loop sites of the VP1 protein [27], was highly conserved among most CVA16 genotype B1 strains, except the (KC117317-GD-2009) strain. No amino acid changes were observed in the 52, 103, 251 and 295 positions in the VP1 protein, which were correlated with virulence among all Chinese CVA16 strains belonging to the B1 subgenotype, except for the (EU812514-AH-2008 and KX595295-GD-2014) strains. Mutation at the 161 position in VP2, reported to increase the virulence in mice, was observed in only two CVA16 strains with amino acid mutations (T161Q/S) at this point. Less pathogenic strains dominated, and highly virulent strains were not endemic in China. Mutations at a single point were the most commonly observed in the P1 gene. Large mutations or deletions were not observed in these regions, indicating a low probability of genomic recombination resulting in the production of a new virus. The function and significance of amino acid mutations and other amino acid changes on the VP1–VP4 protein should be closely monitored and requires further investigation.

This study obtained all available P1 coding region sequences from the GenBank database to identify the genetic characteristics of CVA16 strains related to HFMD in China. The B1a and B1b subtypes were the most prevalent strains in China, particularly in coastal areas, and some CVA16 strains might show inter- and intra-typic recombination in the P1 region. The amino acid sequences of the VP1–VP4 protein encoded by the P1 gene were highly conserved, although several amino acid mutations were observed in the structural protein encoded by the P1 coding region of the Chinese CVA16 strains. No vaccine or other treatment is available at present to prevent or treat CVA16 infection. These findings provide a theoretical foundation for clinical prevention and treatment of HFMD caused by CVA16, as well as the basis for future vaccine development. Long-term and continuous surveillance of CVA16-related HFMD is necessary.