ABSTRACT

The molecular epidemiology of CVA16 in China between 1999 and 2008 reflects a pattern of endemic cocirculation of clusters B1a and B1b within subgenotype B1 viruses. The annual evolution rate of CVA16 was estimated as approximately 0.91 × 10−2 substitutions per synonymous nucleotide/year and is slightly lower than that of HEV71.

Coxsackievirus A16 (CVA16) and human enterovirus 71 (HEV71) are the two major causative agents of hand, foot, and mouth disease (HFMD) (12, 17). In recent years, numerous large outbreaks of HEV71-associated HFMD, which were often accompanied by severe complications, including death, increased research interest in HEV71 strains (1, 14, 16). In contrast, little interest was paid to CVA16 strains because CVA16-associated HFMD was usually mild and benign (2, 12). Although CVA16 is genetically most closely related to HEV71, the genetic diversity and molecular evolution of CVA16, unlike those of HEV71, have not been fully described (5, 7, 12, 13). Cocirculation of CVA16 and HEV71 has been proven to have contributed to the serious outbreaks of HFMD that have occurred in China since 2007 (17); therefore, the genetic variability and the evolution of CVA16 were determined in this study.

The 42 CVA16 strains evaluated in this study were isolated from HFMD patients from different geographical locations in the Shandong, Gansu, Inner Mongolia, and Qinghai provinces of China between 2007 and 2008 (see supplemental data). To investigate the molecular epidemiology of CVA16 in mainland China, 24 additional Chinese CVA16 sequences found in Beijing and Guangdong provinces between 1999 and 2005 and 35 international CVA16 sequences (obtained from the GenBank database) were also analyzed.

A total of 66 Chinese CVA16 sequences were divided into three lineages on the basis of phylogenetic analysis (Fig. 1). A 6.5 to 8.1% nucleotide divergence was found among these three lineages, suggesting that the CVA16-associated HFMD outbreaks in China were a result of the coincident circulation of three genetically distinct viruses.

To determine the molecular epidemiology of Chinese CVA16 strains associated with HFMD epidemics, a phylogenetic dendrogram was constructed with 21 Chinese CVA16 sequences (randomly selected on the basis of their genetic relationships) that circulated during the period 1999-2008 in addition to the 35 international CVA16 sequences that represented two known genotypes (A and B) (13) (Fig. 2).

As in a previous study (13), all CVA16 strains could be grouped into genotypes A and B. The prototype G-10 strain differed from the other strains by 27.5 to 30.2% and clustered separately from all other CVA16 strains, including Chinese CVA16 strains, which clearly belonged to genotype B. This finding was based on the fact that the genetic variation between all other CVA16 strains was less than 13.5%. The sequences in genotype B could be further divided into B1 and B2 subgenotypes with a bootstrap support of 100% (Fig. 2). Chinese CVA16 strains isolated between 1999 and 2008 and the majority of international CVA16 strains isolated between 1997 and 2007 formed subgenotype B1, and the 9 CVA16 strains isolated from Japan and Malaysia between 1981 and 2000 formed subgenotype B2. The nucleotide divergence between subgenotypes B1 and B2 was 11.8%.

Phylogenetic classification based on the complete VP1 region (891 bp) of HEV71 has proved to be useful in tracking genotypes of HEV71-associated HFMD over different temporal and geographical outbreaks (1, 3, 8, 14, 17). Previous studies of CVA16 genetic diversity by Li et al. (12) and Iwai et al. (7) showed that CVA16 strains could be divided into three different clusters, called A, B, and C. However, clusters B and C in their studies correspond to subgenotypes B2 and B1, respectively, in our study, in which a difference of at least 15% in the complete VP1 region of CVA16 strains was used to distinguish genotypes (13). Thus, their B and C clusters should be combined into one genotype.

Subgenotype B1 could be further divided into clusters B1a, B1b, and, possibly, B1c. The nucleotide divergence between clusters B1a and B1b was 6.5%. All Chinese strains isolated between 1999 and 2008 belonged to clusters B1a and B1b, which could also be found in Taiwan, Malaysia, Thailand, Australia, Vietnam, and Saudi Arabia during the same period. This indicated that Chinese strains coevolved and cocirculated with those from neighboring countries. Only one strain from Malaysia (SB16087/MAL/2005) belonged to the probable subgenotype B1c, and the nucleotide divergences between clusters B1c and B1a and clusters B1c and B1b were 7.0% and 6.8%, respectively. Close relatives of this strain have not been found among the known CVA16 strains. This probably reflects the lack of additional strains from Malaysia and its surrounding countries and may represent a surveillance gap in these countries.

Because of the lack of a true “ancestor” strain and the apparent presence of multiple lineages, CVA16 sequences were selected on the basis of their genetic relationships in order to estimate the evolutionary rate (substitutions per nucleotide per year [Fig. 2]) (1, 9). Using the oldest strain in each subgenotype/cluster as a reference, the evolutionary distances were calculated by pairwise comparison according to the Kimura two-parameter method of the MEGA program (11, 15). Both Ks (synonymous substitutions per synonymous site) and Kt (all substitutions per site) values were analyzed, and the evolutionary rate was calculated by linear regression of the genetic distance from the oldest strain versus the year of isolation.

The approximate evolution rate of CVA16 was estimated from the differences in VP1/capsid sequences and the date of specimen collection among strains within cluster B1a and strains within subgenotype B2 (Fig. 2 and Table 1). The average evolution rates for the cluster B1a and the subgenotype B2 were 0.91 × 10−2 (Ks estimate) and 2.49 × 10−3 (Kt estimate), respectively (Table 1), which were slightly lower than the rates calculated for HEV71 (1.35 × 10−2) (1) and poliovirus type 1 (3.36 × 10−2), which is similar to HEV71 (10), although CVA16 strains are broadly distributed geographically.

The evolution of HEV71 and that of CVA16 were similar in the sense that their prototype strains were the sole members of genotype A; eventually, they all gave way to modern strains due to molecular evolution. HEV71 was first identified in 1970 in the United States, after which it appeared as two cocirculated genotypes (genotype B, subgenotypes B1 to B5; genotype C, subgenotypes C1 to C5) during a 40-year course of evolution (1, 6, 14, 17). In contrast, the evolutionary rate of CVA16 is relatively slow in light of the fact that CVA16 was first identified in 1951 in South Africa and that all CVA16 strains, except the prototype strain, formed a single genotype (genotype B) containing two subgenotypes (B1 to B2) as a result of approximately 60 years of evolution.

In conclusion, our study revealed that the molecular epidemiology of CVA16 in China during the years 1999-2008 reflects a pattern of endemic circulation of clusters B1a and B1b within subgenotype B1 viruses; however, in comparison with HEV71, the evolutionary rate of CVA16 is relatively lower.

Phylogenetic dendrogram constructed by using the maximum-likelihood method implemented in the PHYLO_WIN program, version 2.0 (4), based on the alignment of the complete VP1 gene sequences of 66 CVA16 strains (from HFMD patients in the Shandong, Gansu, Inner Mongolia, Qinghai, Beijing, and Guangdong provinces of China) and 1 CVA16 prototype strain (G-10). The bootstrap values from 1,000 pseudoreplicates for major lineages within the tree are shown as percentages. Twelve Chinese CVA16 strains, which were isolated from 12 HFMD cases during 2007-2008 and are indicated by filled circles (•), and 9 strains, which were isolated from 9 HFMD cases in two provinces during 1999-2005 and are indicated by a triangle (▴), were selected for further molecular epidemiological analysis. Abbreviations of provinces of China: GS, Gansu; NM, Inner Mongolia; SD, Shandong; QH, Qinghai; BJ, Beijing; GD, Guangdong.

Phylogenetic dendrogram constructed by using the maximum-likelihood method implemented in the PHYLO_WIN program, version 2.0 (4), based on the alignment of the complete VP1 gene sequences of 21 representative Chinese CVA16 strains and other international CVA16 strains of known genotypes listed in the supplementary material. The bootstrap values from 1,000 pseudoreplicates for major lineages within the tree are shown as percentages. The prototype HEV71 strain (BrCr) was used as an outgroup. Strains indicated by filled circles (•) are CVA16 strains isolated from HFMD cases during 2007-2008. Strains indicated by a triangle (▴) are CVA16 strains isolated from HFMD cases during 1999-2005. Bootstrap values greater than 80% were considered to be statistically significant for grouping. A difference of at least 15% in the complete VP1 region of CVA16 strains was used to distinguish genotypes (13). Abbreviations of provinces of China: GS, Gansu; NM, Inner Mongolia; SD, Shandong; QH, Qinghai; BJ, Beijing; GD, Guangdong.

Estimation of the nucleotide substitution rate in the VP1 capsid region of CVA16

ACKNOWLEDGMENTS

This study was supported by the National Natural Science Foundation of China (project no. 30900063) and the National Key Technology R&D Program of China (project no. 2008BAI56B00, 2009ZX1004-201, and 2009ZX1004-202).

We acknowledge all laboratories that isolated the viruses necessary for this study.