Molecular Epidemiology and Evolution of Human Enterovirus Serotype 68 in Thailand, 2006–2011

Background Publications worldwide have reported on the re-occurrence of human enterovirus 68 (EV68), a rarely detected pathogen usually causing respiratory illness. However, epidemiological data regarding this virus in particular on the Asian continent has so far been limited. Methodology/Findings We investigated the epidemiology and genetic variability of EV68 infection among Thai children with respiratory illnesses from 2006–2011 (n = 1810). Semi-nested PCR using primer sets for amplification of the 5′-untranslated region through VP2 was performed for rhino-enterovirus detection. Altogether, 25 cases were confirmed as EV68 infection indicating a prevalence of 1.4% in the entire study population. Interestingly, the majority of samples were children aged >5 years (64%). Also, co-infection with other viruses was found in 28%, while pandemic H1N1 influenza/2009 virus was the most common co-infection. Of EV68-positive patients, 36% required hospitalizations with the common clinical presentations of fever, cough, dyspnea, and wheezing. The present study has shown that EV68 was extremely rare until 2009 (0.9%). An increasing annual prevalence was found in 2010 (1.6%) with the highest detection frequency in 2011 (4.3%). Based on analysis of the VP1 gene, the evolutionary rate of EV68 was estimated at 4.93×10−3 substitutions/site/year. Major bifurcation of the currently circulating EV68 strains occurred 66 years ago (1945.31 with (1925.95–1960.46)95% HPD). Among the current lineages, 3 clusters of EV68 were categorized based on the different molecular signatures in the BC and DE loops of VP1 combined with high posterior probability values. Each cluster has branched off from their common ancestor at least 36 years ago (1975.78 with (1946.13–1984.97)95% HPD). Conclusion Differences in epidemiological characteristic and seasonal profile of EV68 have been found in this study. Results from Bayesian phylogenetic investigations also revealed that EV68 should be recognized as a genetically diverse virus with a substitution rate identical to that of enterovirus 71 genotype B (4.2×10−3 s/s/y).


Introduction
Human enterovirus 68 (EV68) is a sporadically detected viral pathogen associated with respiratory tract illnesses. EV68 belongs to the family Picornaviridae which contains socially and economically important pathogenic viruses such as foot-and-mouth disease virus, human enterovirus, rhinovirus, and hepatitis A virus. It has been classified in the genus Enterovirus species D which currently comprises 2 important human pathogens, EV serotype 70 and serotype 94. It possesses a positive single-stranded RNA genome of ,7.5 kb in length encased by a highly structured icosahedral capsid. The viral genome is composed of a 59-untranslated region (59UTR) followed by a long single open reading frame and terminated by a short 39UTR with genetically encoded poly-A tract. The coding sequence is co-translationally processed by a virus encoded protease to yield 4 structural (viral capsid protein (VP)-1 to VP4) and 7 nonstructural proteins (2A to 2C and 3A to 3D) which are responsible for viral replication, protein processing and also contribute to shutting down the host cell's protein production.
Many methods have been established for EV typing. Normally, traditional methods rely on biological properties of the viruses such as antigen distinction which subdivides EV into 'serotypes'. However, this method is time consuming, expensive, and could not classify some newly identified EVs. Consequently, as alternative and sensitive methods of classification of all EVs, reverse transcription polymerase chain reaction (RT-PCR) and nucleotide sequencing of the entire genome or partial VP4 or VP4/VP2 regions combined with the results from VP1 gene have recently been utilized for characterizing EVs [1][2][3][4]. EV68 was originally isolated in California in 1962 from children with pneumonia and bronchiolitis by using virus isolation and seroneutralization approaches [5]. In contrast to other EVs, EV68 shares common biological properties with human rhinovirus (RV), which is acid sensitive virus and grow efficiently at an optimum temperature of ,33uC [6,7]. The virus is usually found predominantly in respiratory illness patients. Infection by this virus can cause various disease severities ranging from mild respiratory illnesses such as common cold to severe acute lower respiratory tract infections (ALRTI) including pneumonia, wheezing and bronchiolitis [8][9][10][11]. Three fatal cases caused by severe respiratory illnesses and associated with EV68 infection have been reported [10,11]. Recent studies have suggested that asthmatic individuals infected with EV68 have a propensity to develop unstable asthma or an acute attack [9]. Despite the medical burden of the disease, no effective antiviral therapies have been approved for either prevention or treatment of EV68 infection.
Awareness of its precise global distribution, longitudinal epidemiological profiles, its pathogenic role, and evolutionary history of EV68 are also still lacking. Furthermore, since the re-occurrence of EV68 infection in several countries worldwide [8], only 3 epidemiological studies of EV68 have been reported from the Asian continent (2 from Japan during the study periods of 2006-2009 and 2009-2010 and 1 from the Philippines during 2009-2010) [9][10][11].
To address these concerns, we conducted a prospective longitudinal population-based study among Thai children aged 15 years (yrs.) and below with acute respiratory tract illness (ARTI) symptoms during 2006-2011. Semi-nested PCR amplification covering a highly conserved 59UTR up to and including the 59 terminus of the capsid VP2 gene (59UTR/VP2) and nucleotide sequencing were applied to elucidate many aspects of the epidemiology and evolution of the virus. The primary objective was to determine whether EV68 infection was associated with ARTI among a large number of children in Thailand. The secondary objective was to better understand the evolutionary history of EV68 circulating lineages and to determine how evolution of this virus might have been shaped by its population dynamics.

Study Population
We retrospectively analyzed 1822 stored respiratory specimens obtained from 1785 children with ARTI complications who visited several hospitals located in Bangkok, Thailand between February 14, 2006 andNovember 8, 2011. The main characteristics of these populations are summarized in Table 1. Of these eligible specimens, 24 originated from 12 children and collected in the course of the same hospitalization for treatment follow-up. Twenty-five children had multiple hospitalizations; 22 had 2 admissions and 3 had 3 admissions. Therefore, a total of 1810 respiratory specimens were available for testing.
To assess the distribution of respiratory virus infection with regard to the specific children's age and severity of respiratory illness, the subjects were categorized into 4 groups as follows: infant (#2 yrs.), pre-school (.2 to 5 yrs.), primary-school (.5 to 12 yrs.), and secondary school children (.12 to 15 yrs.). Of 1810 cases, 26 cases (1.4%; 14 outpatients and 12 hospitalized patients) could not provide specific information with regards to age, 39.2% were infants (n = 709), 23.2% were pre-school (n = 420), 26.5% were primary-school (n = 480), and 9.7% were secondary school children (n = 175) ( Table 1). The age distribution of enrolled patients was between 1 day (d.) to 15 yrs. Clinical severity of ARTI cases was defined as follows: pediatric patient with ARTI complications not requiring hospitalization and O 2 supplementation and normal to slightly increased respiratory rate was defined as 'mild' case, patient with ARTI symptoms, requiring hospitalization, and either did or did not require supplemental O 2 was defined as 'moderate' case, patient with severe respiratory failure with O 2 saturation ,90%, requiring medical ventilation or O 2 supplementation, and admission to the Intensive Care Unit (ICU) was defined as 'severe' case. Consequently, 32.1% (n = 581), and 0.9% (n = 16) of the entire study population were considered as moderate and severe cases, respectively, while the remaining cases were mild ARTIs (67%; n = 1213).

Detection of EV68 in ARTI Patients
Of children with ARTI complications, 812 were infected by at least 1 respiratory virus accounting for 44.9% of all enrolled cases. Rhino-enterovirus (RV-EV) was detected in 15.9% (288/1810). Among these, 25 were confirmed as EV68 infection accounting for 8.7% of RV-EV positive specimens and indicating a prevalence of 1.4% in the entire study population. To investigate the relationship between EV68 characterized in Thailand (EV68-TH strains) and the previously defined strains, phylogenetic tree of the 59UTR/VP2 region ( Fig 1A) and their sequence identity matrix ( Table 2) were constructed. The 59UTR/VP2 region of EV68-TH strains exhibited 86.3-97.1% nucleotide identities to each other. EV68-TH strains shared 71-84.1% nucleotide identity with the first identified strain 'Fermon' (AY426531) and 67.7-93.7% with those of strains identified from other countries. All EV68 infected cases were confirmed by amplification of the VP1 gene. Despite performing semi-nested PCR and nested-PCR amplifications yielding a 2 nd PCR product of ,850, 450 or 250 nucleotides (nt), only 11/25 EV68-positive specimens could be amplified. Fig 1B shows  Demographic data and clinical characteristics of these 25 EV68positive children are summarized in Table 3.
Results from RV-EV specific PCR amplification and nucleotide sequencing also revealed that 15 specimens had other EV infections ( Fig 1A). Of these, 1 specimen was classified as poliovirus Sabin-2 with similar proportions of poliovirus Sabin-3 (vaccine strains confirmed), coxsackie virus (CA)-A9, echovirus 9, echovirus 6, and EV109 whereas CA-B5 and CA-B2 were identified in 2 samples each. Five samples were positive for CA-A24. The overall frequency of detection was 0.8% (15/1810).

Clinical Manifestations of EV68 Infected Patients
EV68 were detected as a sole pathogen in 16/25 cases (64%). Co-infections by other respiratory viruses were found in 9 cases (36%) including 1 with respiratory syncytial virus (RSV) type A, 3 with influenza B virus (Flu-B), and the remaining 5 patients were co-infected with pandemic H1N1/2009 virus (pH1N1/2009). As shown in Table 4, age specifics of children with ARTI complications were 1% infant (7 cases; 1 outpatient and 6 hospitalized patients), 2% pre-school (2 cases; all were hospitalized cases), 2.1% primary-school (10 cases; 9 outpatients and 1 hospitalized patient), and 3.4% secondary school children (6 cases; all were outpatients). The average age of EV68-associated ARTI cases was statistically higher than that of the non-infected group (7.6+5.0 yrs. vs 5.0+4.4 yrs.; p = 0.04). Without co-infection by other pathogens, the average age of EV68-associated ARTIs was 5.8+4.4 yrs. The majority of EV68 cases were primary school (40% of EV68 cases and 2.1% of primary school cases) and secondary-school children (24% of EV68 cases and 3.4% of secondary-school cases). Of EV68-positive children, 64% were children aged .5 yrs. (p = 0.016). Female-to-male ratio of EV68 infected patients was 1.5:1. As for the contribution of the virus to disease severity, the majority of positive cases had moderate severity (32% (8/25)). One patient presented with severe clinical characteristics and required admission to the ICU. EV68associated ARTI represented 6.3% (1/16) of all severe cases. The NJ trees were constructed from nucleotide alignments of (A) the 59UTR/VP2 (612 nt) and (B) partial VP1 region (257 nt) using MEGA software. The genetic distances were calculated according to the Kimura-parameter model. Bootstrap support values with 1000 pseudo-replicates were plotted at the internal branch nodes. EV species A to D were labeled. Period of the sample collection was plotted for each lineage. EV68 and other EVs identified in this study (TH-strain) are indicated in black and white circles, respectively. Virus sequence names or location of isolations shown with GenBank accession numbers were reference strains. Scale bar indicates number of nucleotide substitutions per site. doi:10.1371/journal.pone.0035190.g001 None of the patients died. During the study period, hospitalized patients with ARTI were more likely to be infected with EV68 compared to outpatients (1.5% (9/597) vs 1.3% (16/1213)) but this difference did not show any statistical significance. Summary of discharge diagnoses of EV68-confirmed cases included viral pneumonia (n = 8), asthma (n = 1), and influenza-like illness (ILI) (n = 16). Clinical manifestations of the 9 hospitalized patients with EV68 infections are summarized in table 5. Supplemental O 2 and  medical ventilation were required in all EV68-hospitalized cases. The most common clinical characteristics of these patients were fever, cough and dyspnea (77.8%), and wheezing (66.7%).
Hospitalization lasted a median of 5 d (ranging from 3-16 d).
Four patients with EV68 infections had no medical history while 3 had underlying disease such as gastro-esophageal reflux, congenital heart disease, and autism and 2 had a history of asthma. One of these had RSV-A co-infection in the nasal aspiration specimen. However, there was no significant difference in disease severity between the case with and without RSV co-infection.

Seasonality of EV68 Infections
During the first 2 years (February 2006-July 2008), there was no EV68 infection in the sample population (n = 383). EV68 cases were detected in June, July, and September (rainy season) with 0.9% annual prevalence (5/584; all were outpatients). Increased incidence of EV68 infection associated with ARTI was found in 2010 and 2011. In 2010, the virus could be found during a 6-month period with 1.6% annual prevalence (10/611; all were outpatients). During 2011, EV68 infections were investigated only during the rainy season with 4.3% annual prevalence (10/232; 1 was outpatient and 9 were hospitalized patients). Of all EV68 infections detected, 8% were identified from ARTI children admitted to the hospital in summer, 72% during the rainy season, and 20% in winter. Among the meteorological factors, relative humidity was the only factor statistically correlated with the overall number of RV-EV infected patients particularly in 2011 (correlation coefficient (r s ) = 0.61, p = 0.04). None of the meteorological factors correlated with the frequency of EV68 detection in Thailand.

Substitution Rate and Evolutionary Timescale of Circulating EV68 Strains
Evolutionary rates for VP1 gene were measured as the number of nucleotide substitutions per site per year (s/s/y). By using the Bayesian Markov chain Monte Carlo (MCMC) method, the overall evolutionary rate of the VP1 gene of the entire EV68 sequence dataset (dataset 2) was estimated at 4.93610 23 s/s/y with a 95% HPD limit of 4.01-5.85610 23 s/s/y. The evolutionary rate was lower (3.2610 23 s/s/y with a 95% HPD of 2.27-4.25610 23 s/s/y) when the values were calculated with dataset 1 including a more homologous sampling timescale of all EV species. According to our analysis of the evolutionary timescale of EV68-VP1 sequences (Fig 3), the major bifurcation of the currently circulating EV68 strains from their first isolated prototype (Fermon) occurred 66 years ago (time 1945.31 with (1925.95-1960.46)95% HPD). Based on molecular signature analysis of the BC and DE surface loops of the viral determinant VP1 protein, 3 clusters of EV68 were categorized with the posterior probability (pp) value of 1 for clusters 1 and 2, and pp of 0.93 for cluster 3 (Fig 4). The date when each cluster branched off from the most common recent ancestor (TMRCA) occurred in the 1970s (1975.78, (1946.13-1984

Selective Pressure on EV68
To evaluate the adaptive molecular evolution of EV68, overall selective pressure operating on the determinant encoding VP1 was examined by estimating the ratio of non-synonymous (d N ) to synonymous (d S ) substitutions (v = d N /d S ) across the lineages on a codon by codon basis. Selective pressure was defined as follows:

Discussion
Our longitudinal study investigated the prevalence and clinical characteristics of EV68 infection among Thai children with ARTI who visited various hospitals in Bangkok during 2006-2011. In recent years, with the use of highly sensitive molecular techniques, sporadically detected EV68 virus has been recognized as a reemerging respiratory pathogen in several countries worldwide since it has been detected with higher frequency in 2008 [8,12]. EV68 is becoming an increasingly important etiological agent as it was the only pathogen identified in respiratory specimens obtained from fatal cases reported in the retrospective studies of the Philippines and Japan [10,11]. In the present study, the 25 episodes of EV68 infection associated with ARTI were investigat- ed by phylodynamic relationship analysis, which was performed by evaluating nucleotide sequences of the region encompassing part of the 59UTR/VP2 and VP1 gene. Our results suggested that children above the age of 5 yrs. seem to be specifically targeted by EV68. Other studies also reported a similar correlation between patient's age distribution and frequency of EV68 infection [8][9][10][11]13,14]. Despite EV68 infections generally affecting younger children, the association between virus infections and respiratory illnesses in adult patients has also been established [12,13,15]. In concurrence with the earlier reports, infection by this virus appears to be an important cause of viral pneumonia in children requiring hospitalization with common clinical presentations of fever, cough, dyspnea, and wheezing [10,12,13]. Among children with a history of asthma who experienced more severe asthma attacks, EV68 was the sole etiological agent detected indicating its role as an aggravating factor to severe persistent respiratory disease. Furthermore, a strong association of asthma exacerbation and the virus infection has also been reported [9]. Historically, EV68 has been regarded as an uncommon pathogen associated with respiratory illness. However, until now, the full spectrum of illness caused by the infection of this virus has remained unclear. The involvement of EV68 in central nervous system disease has been reported by various US research teams. Khetsuriani et al. discovered the virus from an acute flaccid paralysis case in 2005 [14], while the other group found the virus in cerebrospinal fluid samples taken from fatal meningomyeloencephalitis children in 2011 [16]. Further comparative populationbased studies which include other infectious diseases, healthy controls, or adult hospitalized patients should be conducted to elucidate the causal relationship between EV68 infection and disease progression.
In Thailand, no epidemic status of EV68 infection has so far been reported. Our results showed that EV68 probably emerged in Thailand in 2009. Since that time, EV68 infected cases have more frequently been detected during 2010 and, thus far, displayed the highest prevalence in 2011. This finding supported prior observations that the upsurge of EV68 infection has apparently occurred on a global scale since 2010. Distinct seasonal patterns of EV68 infection in ARTI patients were found in Thailand during the study period. Infection rates were higher during the rainy seasons of 2009 and 2011 while no seasonal pattern was apparent in 2010. The variation in EV68-seasonal profile examined in our study was similar to the profiles of respiratory viruses in tropical countries with a higher average temperature throughout the year and less changes in seasonal temperature [17][18][19][20]. Other studies performed in temperate climates observed seasonality of EV68 in the autumn seasons (Japan 2010 [11], France [21], Netherlands [12], and Italy 2008 [22]). These observations suggested that the seasonal profile of EV68 may vary by season, geographic location, and year.
The VP1 protein is the largest and most exposed surface capsid protein containing the majority of motifs important for interaction with neutralizing antibodies and the cellular receptor required for virus entry. The evolutionary rate, divergence timescale and selective pressure influencing EV68 adaptation should be investigated to provide a better understanding of the genetic diversity and evolutionary relationship of EV68. Results obtained in the present study suggested EV68 could be divided into 3 clusters based on their phylogenetic relationship and the unique molecular signatures within the BC and DE surface loops which was supported by high values of posterior probability. Our suggestion also resembled Meijer's report [12]. Our study also inferred the rate of evolutionary change for EV68 as 4.93610 23 s/s/y. This value is in accordance with the evolutionary rates of other picornaviruses such as EV71 genotype B (4.2610 23 s/s/y) [23,24]. In comparison with other picornavirus members, the rate of evolutionary change estimated for the determinant region of EV68 is faster than the mean evolutionary rate of human parechovirus (2.8610 23 s/s/y) [25], EV71 genotype C (3.4610 23 s/s/y) [24] and hepatitis A virus (9.8610 24 s/s/y) [26] whereas it is lower than the rate of poliovirus (3610 22 s/s/y) [27,28]. Consistently, the high evolutionary rate among these viruses relies on several factors including replicase fidelity, highly error prone viral RNA-dependent RNA polymerases resulting in a misincorporation frequency of 1 per 10 3 to 10 4 nucleotides [29], rate of transmission, and any synonymous mutations of viral proteins. These factors increase the number of mutations incorporated in viral genomes overtime and prepare the ground for rapid genetic diversification [30]. Moreover, intra-and inter-serotype recombination is also believed to be a large-impact evolutionary mechanism influencing the genetic and antigenic diversity of enterovirus, poliovirus, and other picornaviruses [31][32][33][34][35][36][37][38][39][40]. Our study also showed that purifying selection plays an important part in shaping the evolution of EV68. This constrained mutation of EV68 can be explained by the limited size and the genetic architecture of the viral genome which is overlapping between structural and functional domains [41][42][43]. With regard to adaptive mutation of EV68 under serological selection, the v values of specific residues in the VP1 were determined. The result suggested that EV68 tends to escape from stabilizing selection by positively mutating some residues in the antigenic epitope, in particular the residue located on the surface BC and DE loops, resulting in discrimination between the 3 distinct phylogenetic lineages. Differences in the classification of EV68 have been reported by Rahamat-Langendoen et al [13]. Two distinct evolutionary lineages comprising old and new clusters were    Despite having performed this study on stored specimens collected since 2006 and utilizing sensitive PCR approaches, some limitations need to be considered. First, it was possible to conclude that ARTI associated with EV68 infection was first introduced to Thailand in 2009 since EV68-infection was not found between 2006 and 2008. Nonetheless, the limited number of specimens towards the end of 2008 and during the first 5 months of 2009 might be an obstacle to our conclusion. Second, although we utilized specimens which we collected from several hospitals located in many different regions, at least 80% of children were Bangkok residents. Therefore, the enrolled children might not be representative of the entire country population. Increased surveillance together with improved laboratory diagnostic techniques will help reveal epidemiology, clinical association, and relationship between the evolutionary rate and the worldwide circulation of this virus. Analyses at the genome level would also be required for a better understanding of the role of selective pressure in this virus evolution.
In conclusion, our study provides additional epidemiological and clinical data on EV68 infection in Thailand indicating worldwide re-emergence of this virus with variations in epidemiological profile. Our data highlighted the potential importance of EV68 as a causative agent of severe respiratory illness which causes hospitalization in children with specific age distribution. We also showed that purifying selection is the predominant evolutionary force acting on the capsid VP1 of EV68 in that some amino acid residues on the surface exposed loops of the viral epitope are positively selected for escape mutation.

Study Population and Sample Collection
Posterior oropharyngeal and nasal swabs specimens were collected from non-hospitalized patients who had been diagnosed with ILI. Nasopharyngeal aspirations were collected mainly from immunocompromised patients with ALRTI complications and required hospitalizations. A case of ILI was defined according to PAHO/CDC guidelines for influenza surveillance as a progression of fever (.38uC) and either a symptom of cough, sore throat or pharyngitis. Inclusion and exclusion criteria for ALRTI patient enrolment have been described previously [20]. All samples were collected in viral transport media with the addition of antibiotics (2610 6 U/L of Penicillin G and 200 mg/l of Streptomycin) and transported within 48 hours to the Center of Excellence in Clinical Virology, Faculty of Medicine, Chulalongkorn University, for routine respiratory virus diagnostic testing. All respiratory specimens were divided into aliquots and stored at 270uC until further tested. Children admitted multiple times with more than one month between visits were considered as separate illness episodes.

Ethical Consideration
This study was conducted on specimens collected upon conclusion of routine examinations and stored as anonymous. Permission had been granted by the Director of Chulalongkorn King Memorial hospital. Patient identifiers including personal information (name, address) and hospitalization number were removed from these samples to protect patient confidentiality and neither did they appear in any part of document in this study. The research protocol was approved by the Institutional Review Board (IRB number 329/54), Faculty of Medicine, Chulalongkorn University. IRB waived the need for consent because the samples were de-identified.
Laboratory Diagnosis for Virus-associated ARTI Viral nucleic acid was extracted from clinical specimens using a 96-well viral nucleic acid extraction kit (RBC Bioscience, Taiwan). cDNA synthesis was achieved with the M-MLV reverse-transcription system (Promega, US) and random hexamer primers according to the manufacturer's recommendation. Each specimen was tested for common respiratory viruses including RSV-A and RSV-B, Flu-A (pH1N1/2009, seasonal H3 and H1) and Flu-B [44]. For RV-EV screening, semi-nested PCR using a primer set covering the 59UTR/VP2 region was performed as described elsewhere [20]. The final PCR product of the specimens in which RV implicated was ,540 nt in length. According to the close relationship between RV and EV genome sequences, samples that contained EV could be identified ,650 nt. The PCR products were purified using the PCRExtract&GelExtract Mini Kits (5PRIME, Germany) and sequenced bi-directionally by an automated sequencer (First BASE Laboratories, Malaysia). EV68 implicated specimens were confirmed by performing PCR amplifications of the VP1 gene as previously described [13].

Sequence Analysis
Sequence data for each clinical strain was formatted and assembled by Seqman program of DNASTAR Software (v5.0). To study the distribution and diversity of EV68 identified in this study compared to other EVs, nucleotide sequences of the 59UTR/VP2 were multiple aligned by using ClustalW implemented in the Bioedit program (v7.0.9). Phylogenetic trees were constructed by neighbor-joining (NJ) method implemented in the MEGA program (v5). The reliability of the NJ tree was estimated using 1000 bootstrap pseudo-replications. Phylogenetic distance among sequences was measured using Kimura's two-parameter model.

Meteorological Data
The climate of Thailand is tropical savanna with 3 different seasons including winter (mid-November to mid-February), summer (mid-February to mid-May), and rainy (mid-May to mid-October). Bangkok meteorological data during the study period were provided by the Thailand Meteorology Department (http://www.tmd.go.th/en/). Meteorological data including temperature (degree Celsius (uC)), rainfall (millimeter), and relative humidity (%) were routinely measured at 3-hour intervals in the Bangkok metropolis (standard code 455201) at latitude 13.43.35uN and longitude 100.33.36uE. To determine whether the meteorological factors have an impact on or are associated with EV68 emergences in Thailand, bivariate analysis using 2-tailed Spearman's coefficient (r s ) were conducted using SPSS software (v17.0) (Chicago, USA).

Sequence Dataset Preparation
Datasets of the VP1 gene were constructed to perform the Bayesian phylogenetic investigations for estimating the divergence time and determining TMRCA of the recently circulating EV68 lineages. Nucleotide datasets of the VP1 gene were established as follows: dataset 1 was formed by including almost all available EV68-VP1 sequences retrieved from GenBank database by October 31, 2011 and additional nucleotide sequences identified in this study (EV68-TH strains). The homologous sequences of other EVs comprising EV94 for species D (n = 2), CA-A10, EV90 and EV76 as a representative of species A (n = 4) were included in this data set. EV species B and C were assigned as outgroup. Each isolate was labeled with its corresponding year of sampling or isolating. In addition, an alignment of data set 1 containing 122 sequences corresponding to the VP1 gene was constructed (length = 723 nt). Dataset 2 included only EV68-VP1 sequences. The total sequence number of this dataset was 106 sequences with the same length as data set 1.

Estimated Rate of Evolution and Divergence Time of the Recent Spread of EV68
The rate of substitutions and divergence time of EV68 were estimated by using a Bayesian MCMC approach as implemented in the Bayesian Evolutionary Sampling Tree (BEAST) package (v1.6.2) [45]. The best fit substitution model was chosen by performing a maximum likelihood analysis using the Modeltest package (v3.7) [46]. The best fit model suitable for calculation of nucleotide substitutions was General Time Reversible (GTR)+I (proportion of Invariant site)+G (c-distribution) model, allowing for nucleotide rates to vary among sites within the capsid coding sequence alignments. Alignments were performed by keeping the first and second codon positions in one partition and the third position in a separate partition as (1+2)+3 codon partitions. In order to accommodate rate variation at the third codon position, the relative rate parameter was separately estimated in each partition.
Relaxed lognormal molecular clocks were employed and followed by allowing substitution rate variations among branches on the trees. Molecular model specification was selected using Bayes Factors and used to quantitatively estimate the growth rate and demographic parameter. The dynamic among study populations was estimated by performing a Bayesian Skyline plot. The Bayesian MCMC chain lengths were 10 million generations with sampling every 10000 generations and discarding 10% of the chain as burn-in. Convergence of the chains was achieved by computational run over a sufficient time with inspection of the MCMC samples using TRACER (v1.4). The resulting tree of each run was summarized using Tree Annotator and the maximum clade credibility tree was visualized with FigTree software (v1.1.2).

Measurement of Selective Pressure on EV68
The value of v and the individual site specific selection pressure were measured by using the likelihood based SLAC and FEL contained in the HYPHY package. The database was accessed on the website of Datamonkey interface (http://www.datamonkey. org) [47]. The overall v value was estimated based on NJ trees under the TrN93 substitution model. The significance level for a positively selected site of SLAC and FEL analyses was accepted at ,0.1 (two-tailed binominal distribution). The relative residue abundance within the BC and DE-surface exposed loops were depicted by using WebLogo [48].

Statistical Analysis
Statistical data comparisons between various factors were analyzed by means of Pearson x 2 , unpaired T-test, or Fisher's exact test as appropriate, using SPSS software. All data were considered statistically significant at a p-value below 0.05.