Reemergence of Enterovirus 71 in 2008 in Taiwan: Dynamics of Genetic and Antigenic Evolution from 1998 to 2008

Virus isolation and identification.EV71 isolated in Taiwan in 1986 and from 1998 to 2008 was investigated. The isolates from 1998 to 2008 were isolated from patients at National Cheng Kung University Medical Center in southern Taiwan. Six EV71 isolates from 1986 were isolated by the Virology Laboratory of Kaohsiung Medical University Hospital. Specimens from suspected EV-infected patients were inoculated into appropriate tissue cultures including A549, RD, GMK, and MRC-5 cells. EV71 strains were identified and typed antigenically by either neutralization tests or immunofluorescence tests using monoclonal antibodies (Chemicon International Inc.) (53).

Sequencing of EV71.One-tenth of the EV71 isolates were selected for sequencing analyses each month. The EV71 isolates were from a random sampling of patients with diverse clinical presentations, ranging from uncomplicated hand-foot-and-mouth disease to encephalitis and death. Viral genomic RNA was extracted from cell culture isolates by using a Viral RNA Purification Kit II (Geneaid), followed by reverse transcription-PCR and VP1 sequencing as previously described (7). Full-length genome sequences were determined as follows. Sequencing on both the 5′ and the 3′ termini of the genome was performed by 5′ RACE (RACE) and 3′ RACE Systems (Invitrogen) for rapid amplification of cDNA ends, according to manufacturer's instructions, with primers EV71-86 5GSP1, EV71-86 5GSP2, EV71-86 5GSP3, EV71-86 3GSP1, and EV71-86 3GSP2 (see Table S3 in the supplemental material). The amplified products were cloned into pGEM-T Easy (Promega) and sequenced with T7 and SP6 primers. The EV71 full-length genome product was amplified as previously described (5). Briefly, the viral genomic RNA was extracted, and reverse transcription-PCR was performed by using SuperScript II reverse transcriptase (Invitrogen) for reverse transcription and Advantage 2 polymerase (Clontech) for PCR. PCR products were cloned using a TOPO XL PCR cloning kit (Invitrogen) and sequenced. The specific primers used for genome sequencing are indicated in Table S3 in the supplemental material (18). The sequences were assembled with the ContigExpress module of Vector NTI Advance 8 (Invitrogen). Multiple sequence alignments were performed using Clustal X, version 1.83 (50). Some of the EV71 sequences from 1998 to 2005 were from previous studies (18, 52, 53). Accession numbers are given in Table S1 in the supplemental material.

Phylogenetic analyses and estimation of the dN/dS ratio in VP1 sequences.A phylogenetic tree was estimated by the general time reversible model of PAUP*, version 4.0b (35, 54). Statistical robustness of the 1,000 data sets was analyzed, and the estimation of the significance of branch lengths was determined by the maximum-likelihood method. Single-likelihood ancestor counting (SLAC) and fixed-effects likelihood (FEL) methods at the Datamonkey website (http://www.datamonkey.org ) were performed to examine nonsynonymous and synonymous substitution rates (dN and dS, respectively) for different genotypes and the selection signature in the EV71 VP1 protein coding region (33, 34). Selection pressure by the dN/dS ratio for each VP1 codon was measured, and P values were also calculated for these residues. The cutoff P value (<0.1) for the two-tailed extended binominal test was used to classify a site as positively or negatively selected in the SLAC method. In addition, the cutoff P value (<0.1) for the single degree of freedom likelihood ratio test (a chi-squared asymptotic is used) was used to classify a site as positively or negatively selected in the FEL method. The three-dimensional structure of the EV71 VP1 protein was predicted by using the (PS)² program (Protein Structure Prediction Server; http://ps2.life.nctu.edu.tw/ ) at the Genomic Medicine and Biotechnology Development bioinformatics website (Bioinformatics Core for Genomic Medicine and Biotechnology Development; http://www.tbi.org.tw/about/index.htm ) (12).

Recombination analyses.To analyze EV71 and HEV-A genomes, we used a transition/transversion rate of 10 and a 50% consensus to exclude the poorly conserved sites (10). Resulting alignments were analyzed using bootscan analysis in SimPlot, version 3.5.1, with a neighbor-joining tree algorithm and maximum-likelihood distance model consisting of 1,000 pseudoreplicates (10, 27).

Neutralization test and antigenic cartography.Neutralization tests were performed using antiserum from patients infected by various genotypes of EV71, and samples were assayed in a microneutralization assay with RD cells (29). Titers were determined to 0.5 of a twofold dilution. The tabular neutralization data were analyzed manually and also using antigenic cartography (13, 39, 40, 47). Briefly, antigenic cartography is a way to visualize and increase the resolution of binding assay data, such as microneutralization data. In an antigenic map, the distance between a serum point S and antigen point A corresponds to the difference between the log₂ of the maximum titer observed for serum S against any antigen and the log₂ of the titer for serum S and antigen A. Thus, each titer in a neutralization assay table can be thought of as specifying a target distance for the points in an antigenic map. Modified multidimensional scaling methods are used to arrange the antigen and serum points in an antigenic map to best satisfy the target distances specified by the neutralization data. The result is a map in which the distance between points represents antigenic distance as measured by the binding assay.

Phylogenetic and epidemiologic analyses of the EV71 VP1 protein coding region.In 2008, a reemergent EV outbreak occurred in Taiwan which resulted in 373 severe cases and numerous fatalities (Taiwan CDC database; http://www.cdc.gov.tw ). The total number of EV71 cases in the National Cheng Kung University Medical Center in the same year also dramatically increased to 367 cases (Table 1). To compare the diversity of EV71 isolates in 2008 with those of other outbreaks, phylogenetic analyses were performed of the VP1 sequences of EV71 strains isolated in Taiwan in 1986 and from 1998 to 2008 and of other EV71 sequences derived from the GenBank (Fig. 1). Phylogenetic analyses revealed that genotypes B1, C2, B4, and C4 were the major genotypes circulating during the outbreaks of 1986, 1998, 2000 to 2001, and 2004 to 2005 in Taiwan, respectively. In 2008, genotype B5 became the dominant genotype, indicating that the major genotypes of EV71 continued to change from 1986 to 2008, alternating between B and C. Genotype B4 (isolate N5101-TW98, where N5101 is the isolate number, TW indicates Taiwan, and 98 indicates 1998) isolates were first detected in 1998 and subsequently became the dominant genotype in the outbreak of 2000 to 2001. Early detection followed by subsequent outbreaks also occurred from 2004 to 2005 and in 2008. Genotype C4 (N3340-TW02) and B5 (N2776-TW03 and N2838-TW03) were first isolated in Taiwan in 2002 and 2003, respectively, and became emergent strains in the outbreaks of 2004 to 2005 and 2008, respectively. Therefore, the epidemiological and phylogenetic results showed that EV71 strains which became dominant genotypes in outbreaks were usually found in circulation prior to the occurrence of the outbreak.

• Open in new tab
• Download powerpoint

FIG. 1.

Phylogenetic analysis of the VP1 protein coding region of EV71. Phylogenetic analysis of the VP1 protein coding regions of EV71 isolates from 1986 and from 1998 to 2008 from Taiwan and from reference strains from the GenBank based on the VP1 protein coding region (nucleotides 2439 to 3278; 840 bp). The phylogenetic tree was estimated under the general time reversible model with the program PAUP*, version 4.0b. Bootstrap values (percentage of 1,000 pseudoreplicate data sets) of over 75% supporting each cluster are shown at the nodes. CA16 was included as an outgroup. AUS, Australia; SAR, Sarawak, Malaysia; KOR, Korea; SIN, Singapore.

Antigenic cartography of EV71 by neutralization.Phylogenetic analysis of the VP1 protein coding region showed that EV71 isolates in these outbreaks were clustered into five genotypes, B1 in 1986, C2 in 1998, B4 in 2000 to 2001, C4 in 2004 to 2005, and B5 in 2008. Similar to the immunodominant proteins of EVs (37), genetic evolution of the VP1 protein coding region may contribute to antigenic diversity of the virus. To characterize the antigenic properties of EV71 that circulated in Taiwan, we selected strains representing various genotypes from outbreaks in 1986 and from 1998 to 2008 for microneutralization tests. In our neutralization table, human antiserum from EV71-infected patients during the period of 1998 to 2008 was titrated against isolates clustered in various genotypes (see Table S2 in the supplemental material). These sera obtained from infected individuals were divided into five categories (anti-C2, -B4, -C4, -C5, and -B5) based on the genotype of the infecting virus. A serum sample (sample YFW) from a healthy worker at the National Health Research Institutes was also included. However, the neutralization test table did not clearly show whether the antigenic differences of these EV71 viruses were contributed by different genotypes, and it was difficult to quantitatively analyze and interpret the results of the antigenic evolution of the EV71 isolates. To resolve this issue, antigenic cartography was used to construct an antigenic map from the microneutralization data. In the antigenic map (Fig. 2), the genotype B1 and B4 viruses clustered together while genotype C2 from 1998 was found in another antigenic cluster distinct from genotype B viruses. In comparison to genotype B, genotype C4 viruses were more closely related antigenically to genotype C2. Interestingly, the reemergent genotype B5 viruses in 2003 and 2008 mapped in a separate cluster that was not closer to either genotype B1/B4 or genotype C2/C4 in the antigenic map, indicating that the reemergent EV71 viruses in 2003 and 2008 were antigenically different from the other genotype B and genotype C viruses in the study.

• Open in new tab
• Download powerpoint

FIG. 2.

Antigenic map of EV71 isolates from 1998 to 2008 isolated in Taiwan. The relative positions of strains (colored shapes) and antisera (uncolored shapes) were adjusted such that the distances between strains and antisera in the map represent the corresponding neutralization assay measurements (see Table S2 in the supplemental material) with the least error. The periphery of each shape denotes a 0.5-unit increase in the total error; thus, size and shape represent a confidence area in the placement of the strain or antiserum. The vertical and horizontal lines represent antigenic distance, and because only the relative positions of antigens and antisera can be determined, the orientation of the map within these axes is free. The spacing between grid lines is 1 unit of antigenic distance, corresponding to a twofold dilution of antiserum in the neutralization assay. Isolate number and genotype of viruses or serum numbers and genotype of infected viruses are indicated.

To determine the amino acids which may contribute to the antigenic mapping characteristics, we aligned and compared the deduced amino acid sequences of selected viruses based on the sequences of their VP1 protein coding regions. The sequence comparisons (Table 2) showed that Glu43, Thr58, Thr184, and Ser240 were specific signature amino acids for genotype B including subgenotypes B1, B4, and B5 while Arg22 and Asp31 were represented in only genotype C2. Moreover, with two exceptions (N2838-TW03 and N2776-TW98), genotype B5 had an aspartic acid at position 164 similar to genotypes C2 and C4. Glu164 was observed in B4, suggesting that this residue may play an important role in the antigenic properties of EV71.

Phylogenetic analysis of EV71 VP1.Phylogenetic analyses of the nucleotide sequences of EV71 in the VP1 protein coding region indicated that EV71 continued to evolve from 1986 to 2008. The dominant EV71 strain in every major outbreak was antigenically distinct from the previous outbreak strain, including genotype changes from B1 to C2, C2 to B4, B4 to C4, and C4 to B5. To determine the association between these nucleotide and amino acid sequence mutations, mean dN/dS ratios were calculated by using the SLAC analysis method on the Datamonkey website (Table 3). We divided both the EV71 sequences from our study and the sequences published in the GenBank into genotype B (including B1 to B5) and genotype C (including genotypes C1 to C5) based on phylogenetic analyses. In our results, both genotype B and genotype C showed low mean dN/dS ratios (0.098 and 0.044, respectively), indicating that most of the nucleotide substitutions were synonymous. To further identify the mutations involved in EV71 VP1 evolution, SLAC and FEL analyses were performed to examine the dN/dS ratios of the individual sites in VP1 protein coding regions (Table 3). In genotype C viruses, 24 and 89 negative selection sites were found by using SLAC and FEL methods, respectively; however, neither SLAC nor FEL found evidence of positive selection in the VP1 protein coding region. In contrast, one (SLAC) or two (FEL) positively selected sites were identified in genotype B sequences, and 38 (SLAC) and 84 (FEL) negative selection sites were identified. In SLAC analysis results, codon 145 was determined as a positive selection site, and codons 145 and 98 were determined to be positive selection sites by FEL analysis. The two positive sites were located in the BC (codon 98) and DE (codon 145) loops of EV71, according to the three-dimensional structural prediction of the EV71 VP1 protein by the (PS)² program (12). These results showed a high frequency of synonymous mutations of EV71 and demonstrated that genotype B but not genotype C was positively selected at codon 145 and/or codon 98 in the VP1 protein.

Intratypic and intertypic recombination between EV71 strains in Taiwan and other HEV-A viruses.In addition to the analysis of the capsid protein, which contains determinants for cell binding, we also analyzed the untranslated region and nonstructural protein coding regions of EV71 in 1986 and between 1998 and 2008. The complete genome sequences were aligned with other HEV-A sequences to examine the genetic evolution of EV71 isolates. Similarity plot analysis was first performed using a 50% consensus sequence of each virus in 1998 or of individual viruses isolated from 1999 to 2008.

In 1998, genotype C2 caused a large outbreak associated with severe encephalitis. By bootscan analysis, there were two recombination crossover points with a high χ² value, which supported the presence of EV71 subgenotype C and coxsackievirus A8 (CA8) genome sequences within the genome of subgenotype C2 isolates (Fig. 3A). These results showed a possible recombination event between sequences of the CA8 nonstructural protein coding region and the EV71 structural protein coding region. Isolates of subgenotype B4, which were initially identified in 1998 in Taiwan, caused large outbreaks from 2000 to 2001 and continued to circulate in 2003. Genotype B4 of EV71 circulating in Taiwan from 1999 to 2002 (N7008-TW99, as a representative) were phylogenetically close to genotype B3 in the 5′ untranslated region and P1 polyprotein coding regions and similar to genotype B2 in the P3 polyprotein coding region, each with a high bootstrap value (Fig. 3B). These results revealed the possible evolution of genotype B4 from genotypes B3 and B2.

• Open in new tab
• Download powerpoint

• Open in new tab
• Download powerpoint

FIG. 3.

Dynamic analysis of EV71 genetic recombination from 1998 to 2008. Complete genome sequences of EV71 isolates of various genotypes were analyzed by bootscan analysis, including 1998 viruses (A), N7008-TW99 (B), N2838-TW03 (C), and S0584-TW04 (D). This analysis was calculated by SimPlot, version 3.5.1, using the neighbor-joining tree algorithm (Kimura distance model) in a sliding window of 500 bp with a 20-bp step. The EV71 genetic map is shown at the top of each panel.

After the occurrence of genotype B4 for 5 years (1998 to 2003), only some sporadic genotype B4 and B5 viruses were isolated in 2003 in Taiwan. We also analyzed the full genome sequence of one genotype B5 strain (N2838-TW03) from 2003 (Fig. 3C). Both similarity plot and bootscan analysis showed that the genotype B5 genome was more similar to genotype B4 (92.9 to 97.8%) than to other genotypes, especially in the nonstructural region. These results were supported by a high bootstrap value. In 2004 and 2005, a subgenotype shift was observed from genotype B4 to genotype C4, which emerged and became predominant in Taiwan. Bootscan analysis revealed that genotype C4 was similar to EV71 genotype C in the P1 polyprotein coding region and to genotype B from the 2B to 3B protein coding region (with a high bootstrap value of up to 95%) (Fig. 3D). The patterns from similarity plots and bootscan analysis performed with strains from 2004 and 2005 were similar to those reported with N3340-TW02 in a previous study (18). These results suggested that the recombination of the subgenotype C structural protein coding region with the genotype B nonstructural protein coding region resulted in the subgenotype C4 virus responsible for the 2004 and 2005 outbreaks. In 2008, the results showed that genotype B5 exhibited phylogenetic similarity to genotype B4 in the nonstructural region (data not shown), which was identical to N2838-TW03 in 2003.