Comparative Genomic Analysis of Coxsackievirus A6 Strains of Different Clinical Disease Entities

Background Studies regarding coxsackievirus A6 (CVA6) infection were limited. In Taiwan, outbreaks of CVA6 occurred in 2009 and 2010, respectively, but the clinical manifestations were markedly different. We conducted a study to compare the clinical features and genomic sequence between the two years. Methodology/Principal Findings In 2009 and 2010, 205 patients with coxsackievirus A6 (CVA6) infection were treated at Chang Gung Memorial Hospital. Detailed clinical features were obtained from 126 inpatients, 62 in 2009 and 64 in 2010. Between the inpatients in 2009 and 2010, no statistically significant difference was noted in terms of demographics, length of hospital stay and laboratory data. Significantly more patients in 2009 presented with herpangina (82%) while more patients in 2010 presented with hand-foot-mouth disease (HFMD; 67%) and skin rash beyond the typical sites for HFMD. Complete genomic sequences were determined and compared for three isolates from patients with herpangina in 2009 and three isolates from patients with HFMD in 2010. The complete sequences showed that 2009 and 2010 CVA6 isolates were indistinguishable by partial VP1 genes, but there were 5 unique nucleotide changes in 3′ UTR, and 23 out of 2201 (1%) amino acids were different. 2010 viruses underwent the largest number of amino acid changes in 3CD protein, which is the precursor of both 3C protease and 3D polymerase. Conclusions Since 2008 in Finland, outbreaks of HFMD due to CVA6 were noted internationally. CVA6 of different genetic background may cause different clinical manifestations such as herpangina and HFMD.

After EV71 epidemic in 1998, virus surveillance for enterovirus activity has been monitored by the Centers for Diseases Control (CDC) of Taiwan. CVA6 ranked among the top five serotypes in Taiwan between 2001 and 2010. In 2009, CVA6 was the most common endemic enterovirus, and it ranked the 3 rd in 2010 [12]. In our previous study [8] focusing on CVA6 infection from 2004 to 2009, 76.6% of the patients presented as herpangina, and HFMD accounted only for 12.8%. However, in 2010 we found most children infected with CVA6 presented with oral ulcers, various vesicular skin lesions and even onychomadesis, a picture totally different from what we saw previously. These observations prompted us to conduct a study to figure out the differences from clinical manifestations of the patients with CVA6 infection to complete sequence analysis of CVA6 clinical isolates identified between the years 2009 and 2010.

Ethic Statement
The study was approved by the Institutional Review Board of Chang Gung Memorial Hospital. Since the data obtained in this study were collected from the patients, who just received regular medical management, by retrospective medical charts review, a written consent from the patients was waived.

Definitions
Herpangina was defined as well-characterized vesicular enanthem and then ulcers of the fauces and soft palate with presentation of fever, sore throat, and decreased appetite. HFMD also had oral ulcers but chiefly on the buccal mucosa and tongue, accompanied with typical vesicular rashes most commonly on the extensor surfaces of the hands, feet, knees and buttocks. Leukocytosis was defined as WBC count $17.5610 3 /mL. Viral co-infection was defined as virus other than enterovirus was also detected in viral isolation.

Complete Genome Sequencing of Coxsackievirus A6
Three isolates from patients diagnosed as herpangina in 2009, and three isolates from HFMD patients in 2010 were selected for complete sequencing analysis. Viral RNA was extracted from virus culture using the Viral RNA Extraction Miniprep System Kit (Viogene, Sunnyvale, CA, U.S.A.) according to the protocol recommended by the manufacturer. The cDNA synthesis was performed using RT-primer and M-MLV reverse transcriptase (ReverTra Ace; Toyobo, Osaka, Japan). Because large genome difference was observed between prototype and recent field isolates and no other complete genome was available in GenBank at the time of this study, the whole genome sequencing of CVA6 was determined by PCR-based primer walking using the primers listed in Table S1. The sequences obtained in this study were deposited in GenBank under the accession numbers JQ946050-JQ946055.

Sequence Analysis
Multiple sequence alignment was conducted using the ClustalW multiple alignment program within the BioEdit Sequence Alignment Editor package, version 7.0.9.0. [16]. The coding region sequences were translated to amino acid sequences and the percent identities between pairs of sequences were calculated by BioEdit. The unique nucleotide and amino acid substitutions of CVA6 isolates in 2010 were analyzed by comparing with the CVA6 genomes in 2009. The sequence data of the prototype strain Gdula were also used for comparison and analysis. The phylogenetic tree based on VP1 genes was performed using the neighbor-joining method in the MEGA tree program, with a bootstrap value of 1,000 [17].

Statistical Analysis
For categorical variable, we used chi-square test, and the Mann-Whitney U test was used for continuous variables. P values are two-sided, and they are considered significant if #0.05. All analysis were performed with the software SPSS, version 17.0.

Clinical Aspect
Of the 205 patients included, 99 patients were identified in 2009 and 106 patients in 2010. The median age was 2.4 and 2.6 years for children in 2009 and 2010, respectively; 68% and 55% of the patients in 2009 and 2010, respectively, were less than three years of age. The male-to-female ratio was 1.41 and 1.79, respectively. There was no statistically significant difference for age and gender between the two groups. Figure 1 showed the monthly distribution of these children, with a peak in June and August, respectively. Table 1 illustrates the demographics, clinical manifestations and laboratory data of the 62 and 64 hospitalized patients in 2009 and 2010, respectively. Demographics of the patients as well as the length of hospital stay were comparable for both groups. Leukocytosis and elevated serum C-reactive protein (CRP) level .40 mg/L were noted in a substantial proportion of the patients in both groups but there was no statistically significant difference between the two groups.
Nearly all the patients had fever, with a median duration of 2 days and ranging from 1 to 5 days for patients in 2009 and 0 to 7 days in 2010. Although around 90% of the patients in each year had oral ulcers, 90% of the patients in 2009 had oral ulcers limited to posterior pharyngeal wall and soft palate only, while 43% of the patients in 2010 had ulcers (P,0.001) beyond posterior pharyngeal wall and soft palate. Likewise, the presentation of skin rash and the distribution of clinical diagnoses were significantly different between the patients in two groups (P,0.001 for both). 82% of the patients in 2009 had a clinical diagnosis of herpangina, while 67% of the patients in 2010 presented with HFMD. Moreover, 44% of the patients in 2010 had skin rashes beyond the typical sites for HFMD, including trunk, neck, face, and perioral area. All the patients in these two years recovered uneventfully.
There were seven and five inpatients with other viral coinfection in 2009 and 2010, respectively. The difference of clinical manifestations described above was still significant after excluding the patients with co-infection.

Complete Genome Sequencing of Coxsackievirus A6
During the study period, the genome of prototype strain Gdula was the only sequence available for CVA6 complete genome analysis. By comparison with the prototype strain, the recent CVA6 isolates in 2009-10 had the nucleotide sequence identities of approximately 88% in 59 UTR, and 82-83% in P1 region. However, the sequence identities of 2009-10 viruses were less than 80% in both P2 and P3 regions (Table S2). Recombination between enteroviruses of various serotypes may explain the findings. However, when comparing the strains between 2009 and 2010, we found no obvious difference in similarity and bootscan analysis. Therefore, consensus degenerate primers designed from one prototype strain were limited and the complete genome in this study was sequenced largely by primer-walking strategy, especially for P2 and P3 regions.
To investigate genetic basis of 2010 CVA6 isolates for the changing in clinical presentations, the phylogenetic tree was constructed based on partial CVA6 VP1 genes. As shown in Figure 2, one of 2010 CVA6 strain (TW/391/10) was closely related to the strains previously reported in Taiwan [3,4,7,10]. This new lineage was also similar to recent isolates in China, France and Indian [18,19]. However, the phylogenetic results showed that 2009 and 2010 CVA6 isolates were indistinguishable by partial VP1 genes, although CVA6 infections between these two years were easily discernible from clinical presentations.
To identify the genomic sequences of CVA6 associated with distinct phenotype, three complete genomes of CVA6 from each year were compared with each other. The untranslated regions (UTR) present at both end of CVA6 genome are crucial for translation and replication of viral genome. There were 5 unique nucleotide changes in 39 UTR of 2010 viruses, but no consistent nucleotide change was observed in 59 UTR (Fig. 3). Comparison of complete genomes showed that 23 out of 2201 (1%) amino acids consistently differentiated the 2010 from the 2009 viruses ( Table 2). The viral capsid proteins (VP1-VP4) are known to involve in receptor binding and antigenic property, thus having a role in eliciting immune response. The capsid protein sequences of 2010 viruses differed from 2009 virus at only 1 and 2 amino acids in VP2 and VP3 protein, respectively. Among two viral encoded proteases (2A and 3C proteases) that are responsible for the process of viral polyprotein, 2010 viruses had 4 and 3 consistent amino acid changes in 2A and 3C proteases, respectively. The viral RNA-dependent RNA polymerase (RdRp), denoted 3D polymerase, is required for viral genome replication. 3D polymerase had 13 amino acid changes, which was composed mostly of amino acid changes among mature viral proteins. 3CD protein, which is the precursor of both 3C protease and 3D polymerase, has the protease activity of 3C protein domain but no RdRp activity. 2010 viruses underwent the largest number of amino acid changes in 3CD protein. Besides, low nucleotide (88.360.23%) and amino acid  (Table 3). Notably, three clusters of amino acid changes were identified in 2010 viruses. Two clusters harboring two amino acid changes (underlined), 'KGH' and 'IL', were located at positions 101-103 in 2A, and 56-57 in 3C, respectively. The other cluster included 4 amino acid changes (underlined), 'IDKIKK', was located at positions 165-170 in 3D. Overall, the results suggest that unique nucleotide and amino acid sequences of 2010 viruses may be associated with the important, but as yet unknown, virus functions that lead to altered pathogenesis of CVA6.

Discussion
From the present study we found clinical manifestations of the patients infected with CVA6 were different markedly between the year 2009 and 2010. Most patients in 2009 presented with herpangina. While the patients in 2010, oral ulcers were frequently seen in the sites other than soft palate and pharyngeal wall, which are commonly seen in herpangina, and three-fourths of the patients had additional skin rashes over trunk, neck and face, which are not typical sites for HFMD. These findings were also observed in another hospital in Taiwan, 2010 [9]. They reported that 22% of patients with CVA6 infection had eruptions around the perioral area, 30% had rashes over trunk 6 neck and 6.5% had generalized skin eruptions. Onychomadesis 1-2 months following HFMD were also noted. Onychomadesis were reported to be an characteristic feature for CVA6 infection in Finland in 2008 [3,4], and later in Japan in 2011 [10]. However, in the present study we could not find the description of onychomadesis from the charts review and could not know the exact incidence rate of onychomadesis.
In 2009, 82.3% of patients were diagnosed as herpangina; only nearly 10% of the patients were HFMD, which were consistent with our previous study [8]. In contrast, most patients in 2010 were diagnosed as HFMD. Data from Taiwan CDC revealed that CV A2, A4, A5, A6, and A10 were the most common serotypes implicated in herpangina from 2000 to 2005 [20]. Herpangina outbreak related to CVA6 was also noted in Japan in 2005 [21]. In contrast, CVA6 has been found to be an emerging causative agent for HFMD outbreaks in Finland [4] and Singapore [5] since 2008 and also in India in 2009 [18]. After the outbreak of HFMD in Taiwan in 2010, the outbreaks were also from Japan in 2011 [10] and from the United States in 2012 [11].
Generally, patients with HFMD have a temperature of 38uC to 39uC lasting 1 to 2 days [2]. In the present study, almost all the patients in 2010 had fever, and the fever lasted longer (45% of cases $3 days) and higher (two-thirds $39uC). In addition, a substantial proportion of the patients had leukocytosis and elevated serum CRP.40 mg/L. All these findings indicated a higher disease severity caused by the CVA6 strain of 2010.
The VP1 gene of enterovirus genome has been extensively used in phylogenetic analysis due to high degree of diversity among virus serotypes [22]. However, in the present study, the entire genome complexity was greater than anticipated based on partial VP1 sequences, which may misinterpret the phylogenetic relationship. Furthermore, the genome sequence of the prototype strain Gdula isolated in 1949 is obviously different from current strains in Figure 2. Phylogenetic analysis of CVA6 VP1 genes. Partial VP1 sequences (nt 2957 to 3306 according to Gdula numbering) of 6 CVA6 strains obtained from this study and 135 partial VP1 sequences of reference strains derived from GenBank were used to perform phylogenetic analysis. The phylogenetic tree was constructed using the neighbor-joining method with 1,000 bootstrap replications, as implemented in MEGA version 4.  this study (Tables S2 and S3) [23] Therefore, the complete genome sequencing of current CVA6 and other serotype enteroviruses may contribute to the study of enterovirus evolution and genetic variations.
The capsid protein region of 2010 CVA6 possessed 3 unique amino acid substitutions; however, preliminary data from indirect immunofluorescent assay (IFA) suggested that these mutations are not sufficient to produce obvious antigenic changes when compared with 2009 viruses (data not shown). Most mutations reported in this study occurred in non-structural proteins and their functional domains have been defined in recent years. For example, the residues His-20, Asp-38, and Cys-109 comprising catalytic triad of poliovirus 2A protease are essential for inhibition of cap-dependent translation of host mRNA by cleavage of eukaryotic initiation factor 4G (eIF4G) [24,25]. Poliovirus 3C protease contains a His-40, Glu-71, Cys-147 catalytic triad and a conserved KFRDIR sequence for RNA recognition [26,27]. However, none of the mutations in this study were found among these functional sites (Table 2). Following the 'right hand' structure of 3D polymerase, two mutations of 3D polymerase at positions 8 and 33, and one at position 436 were located in the 'index' and 'thumb' sub-domain, respectively. The 'ring' sub-domain contained the clustered mutation 'IDKIKK'. Besides, mutations at positions 308 in motif B and 342 in motif D were found in the 'palm' sub-domain [28,29]. Extensive interactions between various sub-domains are important for virus polymerase function, including NTP binding, RNA binding and elongation [30]. In addition to coding region, 2010 CVA6 genome also contained 5 unique nucleotide changes in 39 UTR, in which RNA secondary and tertiary structures existed is important for RNA replication [31]. Nevertheless, 39 UTR, as well as all structural and nonstructural proteins were known to participate in virus lifecycle, the functional effects of mutations in various regions of viral genome need to be investigated further by mutational analysis.
Some limitations should be noted in this study. First, it is a retrospective study. There may be some loss of chart record about the pattern and sites of skin rash. Inter-observer variation may exist, too. This could explain why the proportion of perioral rashes in the present study was lower (only 10.6% of the cases) than that in other studies (22% in one study from Taiwan [9] and 41% (facial rash) in the recent outbreak in the US [11]). Additionally, records of long-term follow-up were lacking, therefore we could not analyze the characteristic manifestations such as the nail abnormalities and desquamation weeks to months after acute episodes. Finally, there were only three isolates representing the 2009 and 2010 epidemics, respectively. More genetic information of 3'-UTR and RdRP from more isolates are needed to confirm the findings in the present study and to define the correlation between genetic background and clinical manifestations.

Conclusions
Since 2008 in Finland, outbreaks of HFMD due to CVA6 were noted internationally. In addition to herpangina, CVA6 has been an emerging cause of epidemic HFMD other than EV 71 and    Author Contributions