Molecular Evolution of Enterovirus 68 Detected in the Philippines

Background Detection of Enterovirus 68 (EV68) has recently been increased. However, underlying evolutionary mechanism of this increasing trend is not fully understood. Methods Nasopharyngeal swabs were collected from 5,240 patients with acute respiratory infections in the Philippines from June 2009 to December 2011. EV68 was detected by polymerase chain reaction (PCR) targeting for 5′ untranslated region (5′UTR), viral protein 1 (VP1), and VP4/VP2. Phylogenetic trees were generated using the obtained sequences. Results Of the 5,240 tested samples, 12 EV68 positive cases were detected between August and December in 2011 (detection rate, 0.23%). The detection rate was higher among inpatients than outpatients (p<0.0001). Among VP1 sequences detected from 7 patients in 2011, 5 in lineage 2 were diverged from those detected in the Philippines in 2008, however, 2 in lineage 3 were not diverged from strains detected in the Philippines in 2008 but closely associated with strains detected in the United States. Combined with our previous report, EV68 occurrences were observed twice in the Philippines within the last four years. Conclusions EV68 detections might be occurring in cyclic patterns, and viruses might have been maintained in the community while some strains might have been newly introduced.


Introduction
Human enterovirus 68 (EV68) is a member of human enterovirus D (HEV-D), which belongs to the Enterovirus genus, Picornaviridae family. EV68 was first isolated in 1962 from 4 pediatric patients hospitalized with lower respiratory tract infection, in California [1]. Since then, EV68 has been identified only sporadically. In the enterovirus surveillance conducted in the United States from 1970 to 2005, only 26 strains were reported over 35 years [2]. However, the number of reported EV68 cases has remarkably increased in recent years [3]. Between 2008 and 2009, we detected EV68 in 21 pediatric patients who were hospitalized with severe pneumonia in the Philippines [4]. Increased EV68 detection has been reported from different parts of the world, including the United States [5], Europe [6,7], Africa [5] and Asia [8,9]. However, it remains elusive whether this increasing trend was due to a real increase of cases or increased detection owing to improvement in detection methods such as usage of polymerase chain reaction (PCR). Retrospective analysis of stored respiratory specimens by PCR revealed detection of EV68 every year in the Netherlands since 1996 [6] and in Yamagata, Japan since 2005 [8], with smaller number of positive cases before the increased detections in 2010 at each study site. However, it is still unknown whether those viruses detected each year had been maintained in the community or newly introduced.
Majority of reported EV68 cases were associated with acute respiratory infections including a considerable number of severe cases especially among children. A total of 5 fatal cases have been reported so far, including 2 pediatric patients hospitalized with pneumonia in the Philippines [4,10,11]. In one previous study among the patients with acute respiratory infections (ARI) in the Netherlands, detection rate of EV68 was significantly higher in ARI patients with bronchiolitis compared to those without it [6], which suggested that EV68 might be more frequently associated with severe clinical symptoms. In the study, EV68 detection rate was highest among patients aged 50-59 years and it was the lowest among those aged younger than 10 years. In a study conducted in Thailand, EV68 was more likely to be detected among hospitalized pediatric patients compared to children who consulted outpatient clinics [9], which concurred with the report from the Netherlands. However, the etiological significance of EV68 is yet to be determined. Therefore, the current study aimed at investigating the molecular evolution of circulating EV68 viruses in the Philippines.

Sample Collection
We conducted prospective studies for respiratory viruses in Leyte island, the Philippines from June 2009 to December 2011 which include three studies; pediatric pneumonia study at Eastern Visayas Regional Medical Center (EVRMC) since May 2008, adult pneumonia study at EVRMC since May 2010 and influenzalike illness (ILI) study at outpatient clinics of Leyte Provincial Hospital (LPH), Tacloban City Health Catchment Center (TCHCC) and Tanuan Rural Health Unit (TRHU) since January 2010. Nasopharyngeal swabs were collected from patients hospitalized at EVRMC with the diagnosis of severe acute respiratory infection (SARI) as defined by the World Health Organization [12], or those who visited the outpatient clinic of LPH, TCHCC or TRHU with ILI, as defined as fever .38uC accompanied by cough or sore throat. Clinical samples were collected from 1,187 hospitalized pediatric patients from June 2009 to December 2011, 456 hospitalized adult patients from May 2010 to December 2011, and 3,597 outpatients at LPH, TCHC and TRHU from January 2010 to December 2011 (Table 1).
EV68 sequences obtained from the same pediatric pneumonia study between May 2008 and May 2009 were also included in molecular analysis to define their phylogenetic and evolutionary relationship with those detected in this study. Additionally, in order to detect the mixed infections, collected samples were also tested for panels of viruses, including influenza viruses, human metapneumovirus, respiratory syncytial virus, and human rhinoviruses following our previously reported protocol [13].
After obtaining written informed consent, respiratory samples and clinical information were collected from patients or their guardians at the time of hospital admissions or outpatient consultations. The study was approved by the Ethical Committee of Tohoku University Graduate School of Medicine, Japan and the Institutional Review Board of the Research Institute for Tropical Medicine, the Philippines.

RT-PCR and Sequencing
RNA was extracted from clinical specimens using the QIA Viral RNA Mini kit (QIAGEN) according to the manufacturer's instructions. cDNA was synthesized using random primers (Invitrogen) and MMLV Reverse Transcriptase (Invitrogen).

Sequence Analysis
Sequence analysis was conducted by MEGA (version 5) software. Phylogenetic trees were generated by using neighborjoining method, with maximum composite likelihood as a substitution model.
The divergence time of EV68 was estimated by using Bayesian Markov chain Monte Carlo (MCMC) approach implemented in the BEAST package v1.7.1. Lognormal relaxed clock (uncorrected) was used with a tree prior of coalescent Bayesian Skyline, and the Hasegawa, Kishino and Yano 1985 (HKY85) as nucleotide substitution model with gamma as site heterogeneity model. The MCMC chains were run for 10 million iterations, with subsampling at every 10000. A 10% burn-in was discarded and maximum clade credibility tree was summarized using TreeAnovator v1.7.1, with Bayesian posterior probability (BPP) values for each node. Statistical confidence for each parameter estimate was represented by the 95% highest probability density (HPD) intervals around the marginal posterior parameter means. The MCMC process was inspected by using TRACER v1.5. Maximum clade credibility tree was generated by using FigTree v1.3.1.
The strains of previously studied VP1 sequences which were used for generating Bayesian tree were listed in Table S1.

Patient Information
Patient information including demographics, symptoms, radiological findings, O2 saturation level (SpO2), and clinical outcomes were collected using clinical record form. The O2 saturation level was measured at the tip of first fingers or big toes using PalmSAT2500 (Nonin Medical Inc., Plymouth, MN, USA) following manufacturer's instructions The O2 saturation level measured after the initiation of O2 treatment or giving bronchodilator was excluded from the analysis. Information collected from the patients who were positive for EV68 in 2011 was analyzed together with those collected from EV68 positive patients in 2008-09, by using JMP Pro 9.0.2 software (SAS Institute Inc., Cary, NC, USA).

Patients
Among a total of 33 patients (21 patients in 2008-09, and 12 patients in 2011) who were positive for EV68 in respiratory samples, the median age of the patients was 1 year and 8 months ranged from 1 month to 76 years, and 18 of them were male (54.5%).
In 2 adult patients, cough and difficulty in breathing were observed, however no wheezing was found. Rales on auscultation and infiltrates on chest radiographs were found in both patients. SpO2 levels before O2 therapy in those 2 patients were 98% and 95% respectively.
For the EV68 positive patient who was detected in ILI study, chest indrawing and rales were not found in physical examination.
Among 11 hospitalized patients, 9 were discharged and 2 died, including 1 fatal case from the pediatric pneumonia study and another from the adulthood pneumonia study. The fatal adult patient with EV68 had a comorbidity of hepatic cirrhosis.

Phylogenetic Analysis
The viral genome of VP1 region was amplified in 7 out of the 12 samples positive for 59UTR, and VP4/VP2 region was amplified in 4. The samples positive for VP1 and VP4/VP2 were collected from the pediatric pneumonia study. Detected VP1 and VP4/VP2 sequences were deposited in GenBank with accession number of AB817710-AB817716, and AB829882-AB829891 respectively. Phylogenetic trees were drawn based on 59 UTR, VP1, and VP4/ VP2 sequences including those detected between 2008 and 2009 in the Philippines (Figure 2).
Of the sequences obtained from 9 hospitalized children in 2011, 5 (TTa-11-Ph245, 266, 269, 344, and 395) belonged to the same cluster as those from 19 hospitalized children in the Philippines between 2008 and 2009 based on a phylogenetic tree of 59UTR (Figure 2a). Of the remaining 4 sequences, 2 (TTa-11-Ph224 and 272) formed a distinct lineage together with sequences obtained from 2 adult hospitalized patients (TTaAP-11-Ph121 and 129) and 1 ILI patient (ILI-11-Ph227) in the present study. This lineage also included sequences detected from Senegal, England, and New York. The remaining 2 (TTa-11-Ph250 and 257) belonged to the cluster including previously reported sequences from the United States (TX03), the Philippines (TTa-08-Ph561) and South Africa (SA1354) (Figure 2a).
Among 12 patients who were positive for EV68 in respiratory specimens collected in 2011 (9 from the childhood pneumonia study, 2 from the adult pneumonia study, and 1 from the ILI study), VP1 sequences were amplified only for 7 samples from childhood pneumonia study. On the phylogenetic tree based on VP1 sequences, 5 out of 7 sequences detected in 2011 were in the same lineage; lineage 2, along with sequences detected from the Philippines in 2008 (TTa-08-Ph343, 451, 513, 519, 560, 597, and    deletions at the positions 2806-2808 in the VP1 region compared to the Fermon strain (the prototype strain of EV68), which resulted in codon deletions [16]. Although 3 strains from the Philippines shared the same branch, those identified in 2011 were clustered in a different branch from the strain in 2008 (Figure 2b).
One of the strains detected in the Philippines in 2011 (TTa-11-Ph257) belonged to lineage 2 on the phylogenetic tree based on VP1 sequences, however, its 59UTR sequences were closely related to previously reported sequences including TX03, SA1354 and TTa-08-Ph561 (Figure 2a). Sequences of these 3 strains (TX03, SA1354 and TTa-08-Ph561) belonged to lineage 3 and formed one cluster together with those detected in Yamagata, Japan in 2006 and 2008 (Y06-2311 and Y08-1737) on the phylogenetic tree for VP1.
Another strain detected in the Philippines in 2011 (TTa-11-Ph250) also belonged to the same cluster with TTa-08-Ph561 on the phylogenetic tree based on 59UTR sequences, however, its VP1 sequences were not amplified. Strains detected in adult pneumonia study (TTaAP-11-Ph121 and 129) and ILI study (ILI-11-Ph227) formed a branch together with other strains including those identified in childhood pneumonia study in 2011 (TTa-11-Ph224 and 272) on the 59UTR tree, however, VP1 sequences of these 3 strains were also not amplified.

Bayesian Tree Analysis
The Bayesian evolutionary tree was generated, in which VP1 sequences of EV68 detected in recent years formed 3 distinct lineages. Minor differences were observed in the estimated divergence time of each node compared to the previous report [9] which might be due to the differences in the data set used for analysis. Among 7 sequences obtained from pediatric pneumonia patients in the Philippine in 2011, 5 (TTa-11  .80 of 95% HPD with median of 2009.39). In the other clusters of lineage 2 which did not include any sequences detected from the Philippines, consisted of those detected in Yamagata; Japan, England and the Netherlands. All the sequences detected in Yamagata formed a distinct cluster which diverged from other sequences detected in European countries.
Sequences of the remaining 2 cases from the Philippines in 2011 (TTa-11-Ph224 and TTa-11-Ph272) belonged to lineage 3 together with sequences including those detected in the Philippines in 2008 (TTa-08-Ph561). Lineage 3 consisted of sequences detected from different parts of the world including Japan, the Philippines, New Zealand, the United States of America, the Netherlands, England, Gambia, and Senegal, and each cluster in this lineage did not correlate with the specific geographical distribution. TTa-11-Ph224 and TTa-11-Ph272 formed a distinct lineage which diverged from the cluster consisted of sequences detected from New York in 2009. Based on this tree, TTa-08-561 detected in the Philippines in 2008 was not likely to be ancestor to TTa-11-Ph224 and TTa-11-Ph272. In lineage 3, sequences diverged into two clusters in 2000 (2000.40-2002.18 of 95% HPD with median of 2000.41), and one of them was consisted of only 2 sequences; TTa-08-Ph561, and Y08-1737, the latter of which was detected in Yamagata, Japan in 2008.

Discussion
We previously reported the detection of EV68 from 21 pediatric pneumonia patients in Leyte island, the Philippines between 2008 and 2009 [4]. In this study, we detected EV68 at the same study site between June and August 2011. Even though the overall detection rate was low (0.23%), all cases were detected over a period of 3 months. Moreover, the detection rate among hospitalized cases was higher (9/1,187, 0.76% of pediatric cases and 2/456, 0.44% of adult cases) than outpatients (1/3,597, 0.028%) (p#0.0001), suggesting that EV68 was more likely to be associated with severe respiratory infections. The hospitalized patients positive for EV68 presented severe respiratory symptoms, and among pediatric patients, the SPO 2 values were considerably lower than the reference value (95%-100%) of healthy infants and children [17]. Furthermore, unlike the previous study from Thailand, interestingly we did not detect any mixed infection with aforementioned viruses in our clinical samples [9].
Combined with our previous report, increased detections of EV68 were observed twice during the 4 year study period, with an interval of approximately 2 years, and no positive cases were detected between those two periods with EV68 detections (Figure 1). Possibly, EV68 recurrence follows a cyclic pattern as suggested by the retrospective analysis of stored samples in Japan and the Netherlands [6,8]. However, in these studies, a small number of cases were detected even between periods with increased EV68 detections.
In 2008 and 2009, majority of the detected EV68 belonged to lineage 2, and only one strain was classified into lineage 3 on the phylogenetic tree based on VP1 sequences. In 2011, most of the detected EV68 shared high similarities in VP1 sequences and also classified into lineage 2, however, two strains of lineage 3 were also detected. These results indicated that EV68 of lineages 2 and 3 were co-circulating in the study area during both 2008-09 and 2011. In the VP1 amino acid sequences detected from the Philippines, nonsynonymous substitutions were commonly found in BC-and DE-loop regions, and unique amino acid sequences in those regions were associated with the categorization of detected strains into different lineages ( Figure S1), which was in line with the previous report from Thailand [9]. It was previously reported that genotype replacement is associated with occurrence of Enterovirus 71 (EV71) outbreaks [18]. In our study, such replacement of EV68 lineages was not observed between two periods of increased detection, therefore, the mechanism for increased detections is still unclear. Observation of longer period is required to reveal the relationship between circulating lineages and increased detections of EV68.
One of the strains detected in the Philippines in 2011 (TTa-11-Ph257) belonged to lineage 2 on the phylogenetic tree based on VP1 sequences, however, its 59UTR sequences were closely related to strains whose VP1 sequences were classified as lineage 3 (TX03, SA1354 and TTa-08-Ph561). This variation may be due to the recombination which is common among members of Picornaviridae family [19,20]. However, we could not make any conclusions since we only amplified partial genome of 59UTR and VP1 of this strain. The role of recombination in evolution of EV68 is yet to be defined.
Among EV68 detected in the study, 5 in lineage 2 formed a cluster which diverged from the cluster of EV68 detected in 2008-09 in the Philippines, on the Bayesian tree ( Figure 3). However, 2 strains from the Philippines which belonged to lineage 3 were not diverged from viruses detected in the Philippines in 2008 and they were more closely related to viruses found in New York, the United States. This result suggests that these strains might have been newly introduced into the community. For other respiratory viruses, strains migration at global level is well documented. For examples, molecular epidemiology of respiratory syncytial virus (RSV) indicated that specific genotype BA of subgroup B has widely spread to different parts of the world [21,22,23]. Despite the geographical distance between New York and Leyte Island, the Philippines, it is still possible that EV68 transmission had occurred at global level. However, number of analyzed sequences in lineage 3 was too small to draw any conclusion.
In conclusion, despite the low detection rate, our study indicated etiological importance of EV68 as a potential causative agent of severe respiratory infections. In addition, increased detections of EV68 in one area were identified twice with a 2-year interval, which suggests that viruses might have been maintained in the community. Given the fact that most of the detected sequences in 2011 were closely associated with those detected in 2008, cyclic patterns of EV68 detections might be due to the accumulation of susceptible population in the community. Our EV68 in ARI Patients, the Philippines PLOS ONE | www.plosone.org study also indicates that transmission of EV68 at global level and introduction of new strains into the community might be an alternative mechanism for increased detections of EV68. Continuous monitoring is thus required to reveal the etiological importance of EV68 and the mechanisms of its emergence in cyclic patterns in the community. Figure S1 Identities in the amino acid sequences of VP1 detected in the Philippines. The amino acid sequences of VP1 region detected from the Philippines were aligned in reference to the protostrain; Fermon (AF081348). The codon positions where the Philippines sequences had identical amino acid residue to Fermon strain was indicated with dots. The BC-and DE-loop regions were indicated in boxes of dotted lines. (EPS) Table S1 The strains whose VP1 sequences were used for analysis. A total of 171 strains, including 15 from the Philippines were used for molecular analysis. The collection year, location, and the accession numbers were obtained from GenBank. (DOCX)