ABSTRACT
In August 2014, an outbreak of enterovirus D68 (EV-D68) occurred in North America, causing severe respiratory disease in children. Due to a lack of complete genome sequence data, there is only a limited understanding of the molecular evolution and epidemiology of EV-D68 during this outbreak, and it is uncertain whether the differing clinical manifestations of EV-D68 infection are associated with specific viral lineages. We developed a high-throughput complete genome sequencing pipeline for EV-D68 that produced a total of 59 complete genomes from respiratory samples with a 95% success rate, including 57 genomes from Kansas City, MO, collected during the 2014 outbreak. With these data in hand, we performed phylogenetic analyses of complete genome and VP1 capsid protein sequences. Notably, we observed considerable genetic diversity among EV-D68 isolates in Kansas City, manifest as phylogenetically distinct lineages, indicative of multiple introductions of this virus into the city. In addition, we identified an intersubclade recombination event within EV-D68, the first recombinant in this virus reported to date. Finally, we found no significant association between EV-D68 genetic variation, either lineages or individual mutations, and a variety of demographic and clinical variables, suggesting that host factors likely play a major role in determining disease severity. Overall, our study revealed the complex pattern of viral evolution within a single geographic locality during a single outbreak, which has implications for the design of effective intervention and prevention strategies.
IMPORTANCE Until recently, EV-D68 was considered to be an uncommon human pathogen, associated with mild respiratory illness. However, in 2014 EV-D68 was responsible for more than 1,000 disease cases in North America, including severe respiratory illness in children and acute flaccid myelitis, raising concerns about its potential impact on public health. Despite the emergence of EV-D68, a lack of full-length genome sequences means that little is known about the molecular evolution of this virus within a single geographic locality during a single outbreak. Here, we doubled the number of publicly available complete genome sequences of EV-D68 by performing high-throughput next-generation sequencing, characterized the evolutionary history of this outbreak in detail, identified a recombination event, and investigated whether there was any correlation between the demographic and clinical characteristics of the patients and the viral variant that infected them. Overall, these results will help inform the design of intervention strategies for EV-D68.
INTRODUCTION
Enteroviruses (EVs) are members of the family Picornaviridae (order Picornavirales) of RNA viruses, with the genus Enterovirus comprising 12 species, designated A to H and J, as well as rhinovirus A to C (1). All EVs are characterized by a single positive-strand RNA genome of approximately 7,500 nucleotides (nt) in length. To date, five types of Enterovirus D (EV-D) species have been described: EV-D68, together with EV-D70 and EV-D94, is an important human pathogen (2–5), while EV-D111 and EV-D120 were only recently assigned (4, 6, 7). EV-D68 was first isolated from hospitalized children with pneumonia and bronchiolitis in California in 1962 (8), and rhinovirus 87 was later reclassified as EV-D68 based on phylogenetic analysis, cross-neutralization, and acid sensitivity (9–11). The EV-D68 genome includes a single open reading frame which encodes four structural proteins (VP1 to VP4) and seven nonstructural proteins (2A to 2C and 3A to 3D), a 5′ untranslated region (UTR) with a hairpin-loop secondary structure, and a 3′ UTR with a poly(A) tract (1).
EV-D68 was detected sporadically from the 1970s through the early 2000s (12–18). However, an increased number of infections by EV-D68 has been reported worldwide during the past decade (15, 16, 18–33). Almost all documented EV-D68 cases during this time were associated with acute respiratory infections, ranging from mild upper respiratory tract illness and asthma to severe bronchitis and pneumonia, with the exception of four isolated cases associated with neurological syndromes, notably, acute flaccid myelitis (AFM), in which viruses were isolated from the cerebrospinal fluid or nasopharyngeal swabs (12, 31, 34). The 2014 outbreak of EV-D68 was the largest outbreak recorded in North America, with more than 1,153 confirmed EV-D68 cases in 49 U.S. states and the District of Columbia (35). The majority of cases were associated with severe respiratory illness in children (36, 37) leading to hospitalization. In addition, there was a surge in EV-D68 infections associated with AFM and encephalitis (36, 38, 39). EV-D68 infections were reported in Europe during the same time period, with three cases associated with AFM (40, 41).
Previous studies have revealed the presence of three major clades of EV-D68, designated A, B and C, which have circulated or cocirculated during different time periods in different geographic regions (18, 38). Most of the viruses sampled from the 2014 EV-D68 outbreak in the United States have been assigned to a novel B1 subclade, especially those sampled from AFM cases (38). Along with high rates of nucleotide substitution, recombination is an important way in which genetic diversity can be generated in enteroviruses (33). However, no recombination events have yet been reported in EV-D68, in part reflecting a lack of full-length genome sequences for analysis (22).
There is growing interest in understanding the epidemiology of EV-D68 in the United States, particularly given its association with severe disease outcomes during the 2014 outbreak. Indeed, little is known about the molecular evolution of the virus during a single outbreak, nor whether differences in clinical manifestation are associated with specific genetic variants. The most important limitation here has been the lack of complete genome sequences of EV-D68 in the public domain, even though such information will assist the design of effective intervention and prevention strategies and help formulate modalities of future treatments. To help determine the evolution of EV-D68 within a single geographic region during a single outbreak, we performed full-length genome sequencing of viruses sampled from children seen at Children's Mercy Hospital—the first hospital reporting the EV-D68 outbreak in 2014—in Kansas City, MO, during August to September 2014 (36). As well as determining the extent and pattern of viral genetic diversity, we screened these genomic data for recombination and assessed whether there was any association between genetic variation and specific demographic and clinical features of the infected patients.
MATERIALS AND METHODS
Study population and patient data collection.Patients, aged 0 to 17 years, admitted to Children's Mercy Hospital from 1 August to 15 September 2014 with a positive test result for enterovirus/rhinovirus by FilmArray respiratory panel assay (Biofire LLC, Idaho) were retrospectively tested for the presence of EV-D68 by real-time PCR (37). EV-D68-positive patients admitted to the pediatric intensive care unit (ICU) were age matched with two patients not requiring ICU care. Data were retrospectively obtained from chart abstraction and entered into a standardized data collection instrument. This project was approved by the Children's Mercy Hospital Institutional Review Board.
A single-term, previously healthy infant with EV-D68 was enrolled in an infant birth cohort based at Vanderbilt University Medical Center and was included as a preemergent isolate in a region geographically close to Kansas City. This is a population-based birth cohort, and respiratory illness surveillance was performed every 2 weeks during the winter season. All information was prospectively collected, and parents gave their informed consent for study inclusion (42).
RNA extraction and reverse transcription-PCR (RT-PCR).Extraction of viral RNA was performed at the J. Craig Venter Institute (JCVI), Rockville, MD, with 140 μl of respiratory samples (nasal swab or nasal aspirates/wash) in transport medium using a QIAamp viral RNA minikit (Qiagen, Hilden, Germany)/ZR96 viral RNA kit (Zymo, Irvine, CA) hybrid protocol. In brief, specimen lysis was performed in Qiagen buffer AVL in a 96-well deep-well plate. Lysate was transferred to a ZR96 spin plate (Zymo), and samples were processed according to the manufacturer's protocol. The cDNA was generated using SuperScript III reverse transcriptase (Thermo Fisher Scientific, MA) from 4 μl undiluted RNA and either oligo(dT)20 or a 1:1 mix of two primers specific to the 3′-terminal region (D68_7333AR and D68_7333BR) (Table 1).
EV-D68 primers used for three strategies of PCR amplification
Three independent approaches were taken when performing PCR amplification of the viral genome. (i) Full-length genome amplicons containing the 3′ poly(A) tail (FL-A) were generated from the cDNA generated by oligo(dT)20 using D68_1F and M13-dt18 (reverse primer). (ii) Full-length genome amplification excluding the 3′ poly(A) tail (FL) was performed using gene-specific RT mix primers and PCR primer D68_1F and a 1:1 mix of D68_7333AR and D68_7333BR. (iii) Generation of complete viral genomes through the amplification of two overlapping amplicons using gene-specific RT was achieved as follows. A small (S) 904-bp amplicon encompassing the EV-D68 5′ UTR was generated directly from RNA using the Qiagen one-step RT-PCR kit (Qiagen) according to the manufacturer's protocol using primers D68_1F and D68_904R. A large (L) 6.8-kbp amplicon encompassing the rest of the viral genome was generated from primer D68_536F and a 1:1 mix of D68_7333AR and D68_7333BR. Amplicons were verified on 1% agarose gels or via the QIAxcel advanced (Qiagen) capillary gel electrophoresis DNA screening platform, and excess primers and deoxynucleoside triphosphates (dNTPs) were removed by treatment with exonuclease I (New England BioLabs) and shrimp alkaline phosphatase (Affymetrix, Santa Clara, CA, USA) at 37°C for 60 min, followed by incubation at 72°C for 15 min. Amplicons were quantified using a SYBR green double-stranded DNA (dsDNA) detection assay (SYBR green I nucleic acid gel stain; Thermo Fisher Scientific), and all four amplicons per genome were pooled in equal concentrations.
EV-D68 complete genome sequencing assembly and annotation.Illumina libraries were prepared using the Nextera DNA sample preparation kit (Illumina, Inc., San Diego, CA, USA) with half-reaction-mixture volumes as described previously (43, 44). For samples requiring extra coverage, the Ion Torrent Personal Genome Machine (PGM) (Thermo Fisher Scientific) was used in addition to Illumina sequencing, in which 100 ng of pooled DNA amplicons was sheared for 7 min and Ion Torrent-compatible barcoded adapters were ligated to the sheared DNA using the Ion Xpress Plus fragment library kit (Thermo Fisher Scientific) to create 400-bp libraries. Sequencing was performed on the Ion Torrent PGM using 316v2 or 318v2 chips (Thermo Fisher Scientific).
Sequence reads were sorted by barcode, trimmed, and de novo assembled using CLC bio's clc_novo_assemble program (45). The resulting contigs were searched against custom, full-length EV-D68 nucleotide databases to identify the closest reference sequence. All sequence reads were then mapped to the selected reference EV-D68 sequence using CLC bio's clc_ref_assemble_long program (46). At loci where both Illumina and Ion Torrent sequence data agreed on a variation compared with the reference sequence, it was updated to reflect the difference. A final mapping of all next-generation sequences to the updated reference sequences was performed with CLC bio's clc_ref_assemble_long program (46). Curated assemblies were validated and annotated with the viral annotation software Viral Genome ORF Reader (VIGOR) 3.0 (47) before submission to GenBank. VIGOR was used to predict genes, perform alignments, ensure the fidelity of open reading frames, associate nucleotide polymorphisms with amino acid changes, and detect any potential sequencing errors. The annotation was subjected to manual inspection and quality control before submission to GenBank.
Phylogenetic analysis.We analyzed 59 newly acquired complete genome sequences of EV-D68 (57 collected in Kansas City during the 2014 outbreak; one preoutbreak sequence collected in Nashville, TN, USA [Vanderbilt University Medical Center], in 2012; and one prototype Fermon virus strain purchased from ATCC) together with all EV-D68 sequences available on GenBank (http://www.ncbi.nlm.nih.gov/GenBank/). In total, we utilized two global data sets of EV-D68: a 110-complete-genome data set (59 JCVI-sequenced sequences and 51 background sequences) and a 357-complete-VP1-gene data set (59 JCVI-sequenced sequences and 298 background sequences). Sequences within these two data sets were aligned separately using the MUSCLE program in MEGA 6.0 with manual adjustment (48). Maximum likelihood (ML) phylogenies (49) were inferred in MEGA 6.0. A general time reversal (GTR) substitution model with a gamma distribution of among-site rate variation and a proportion of invariable sites (GTR+Γ+I) was selected as the best-fit model by ModelTest in MEGA 6.0 and used in all tree inference methods. In addition, trees were inferred using the Bayesian Markov chain Monte Carlo (BMCMC) method available in MrBayes version 3.2.5 (50), run for 1 × 108 steps. Trees were sampled every 1 × 104 steps, with the first 1,000 trees discarded as burn-in. The robustness of the ML tree was assessed by bootstrap analyses of 1,000 pseudoreplicates and by comparison with the topologies sampled in the Bayesian analysis. All phylogenies were rooted with the oldest EV-D68 sequence in GenBank, the Fermon strain, collected in 1962 in California, USA (2).
Inferring the events of virus introduction to Kansas City.A total of 5,000 trees were evenly subsampled from the posterior distribution of trees produced during the BMCMC analysis described above using the LogCombiner program within the BEAST package (51). Each sequence was assigned a discrete geographic state—either “Kansas City” or “other place”—according to our data or its record in GenBank. A parsimony procedure (52) was then used to infer ancestral geographic states given these data and hence to determine the frequency of state change from the “other place” character state to the “Kansas City” character state in each tree. Such changes in geographic state are indicative of independent viral introductions into Kansas City. The mean and 95% confidence intervals of the frequency of introduction events were summarized from the counts in the 5,000 subsampled trees.
Analysis of EV-D68 recombination.Potential recombination within the VP1 and complete genome sequences of EV-D68 was screened using seven methods (RDP, GENECONV, MaxChi, Bootscan, Chimaera, SiScan, and 3Seq) implemented in the Recombination Detection Program version 4.46 (RDP4) (53). Phylogenetic incongruence between different regions and with P values of less than 1 × 10−4 in all methods was taken to represent strong evidence for recombination. To confirm these putative recombination events, we utilized a smaller data set including the recombinant and the parental strains determined above, employing similarity plots and Bootscan analysis as implemented in Simplot version 3.5.1 (54), with a window size of 300 nucleotides (nt) and a step size of 10 nt. Recombination breakpoints were inferred based on the distribution of informative sites supporting the two incongruent tree topologies that maximized the chi-square (χ2) sum (55). The midpoint of the breakpoint region was used to partition the sequence alignment for separate phylogenetic inferences.
Assessing the association between virus phylogeny and the demographic and clinical features of EV-D68 patients.Demographic and clinical information was available for 57 outbreak samples from Kansas City. To determine whether there was phylogenetic clustering by the age and gender of the patient, we grouped viruses into two age classes (<5 years and ≥5 years old) and by gender separately. Similarly, each patient was classified as either positive or negative for each clinical symptom individually. The clinical symptoms analyzed were (i) presence in pediatric ICU, (ii) medical history of asthma or recurrent wheezing, and (iii) requirement for ventilation. The strength of association between the phenotypic features described above and the EV-D68 phylogeny was determined using two phylogeny-trait association statistics, the parsimony score (PS) and the association index (AI) tests, both of which were implemented in the Bayesian tip association significance testing (BaTS) program (56). A null distribution of these statistics was determined using the posterior distribution of trees obtained from the MrBayes analysis described above.
Nucleotide sequence accession numbers.All sequences generated as part of this study were submitted to GenBank as part of the BioProject identifiers PRJNA266349 and PRJNA270340 and assigned accession numbers KT347223 to KT347280 and KT725431.
RESULTS
High-throughput full-length genome sequencing of EV-D68 from the 2014 outbreak.A total of 62 EV-D68-positive samples were chosen for complete genome sequencing, including 60 clinical samples from the 2014 outbreak collected from Children's Mercy Hospital, Kansas City, MO; one preoutbreak sample from Vanderbilt University Medical Center, Nashville, TN, collected in 2012; and one sample from the American Type Culture Collection (ATCC), VR-1197. Oligo(dT)20 and/or two primers specific to the 3′-terminal region of the virus were used for reverse transcription (Table 1). The general pattern of primer sites and the locations of primer targets in EV-D68 genomes are shown in Fig. 1A. To achieve complete genome sequencing, three approaches were attempted from respiratory samples: (i) full-length genome with poly(A) tail (FL-A), (ii) full-length genome without the poly(A) tail (FL), and (iii) two overlapping amplicons—large (L) and small (S) (Fig. 1A). Representative PCR products for each method are shown in Fig. 1B. PCR products of the expected size and complete genome assemblies were obtained with variable success rates: 53% and 58% for FL-A and FL, respectively, and 95% (59 out of 62) for the two-amplicon method, which is superior to a recently published method based on sequence-independent amplification (57). These PCR products were subjected to next-generation sequencing using both the Illumina Mi-Seq and Ion Torrent sequencing platforms. In total, we obtained 59 complete genome sequences of EV-D68, including 57 outbreak samples from Kansas City, one preoutbreak sample from Nashville, and one historical sample from ATCC.
High-throughput sequencing of the complete genome of EV-D68. (A) Schematic representation of the different approaches used to sequence the complete EV-D68 genome. FL-A, full-length genome with poly(A) tail; FL, full-length genome without poly(A) tail; S, small overlapping amplicon; L, large overlapping amplicon. (B) Three examples of the EV-D68 RT PCR product—2014 outbreak strain MO49/2014, preoutbreak strain Vanderbilt N0051U5/2012, and ATCC strain VR1197—resolved by agarose gel electrophoresis and visualized by ethidium bromide staining. Sizes of molecular markers (in base pairs) are indicated on the left side of the gel.
Our sequencing resulted in an average of 62,256 reads/sample with an average coverage of 1,569× (minimum average coverage of 219× and maximum average coverage of 3,324×). Deep sequence analysis revealed an average of 4.54 single nucleotide polymorphisms (SNPs) per virus that are present in more than 3% of reads and hence above the background level of SNPs expected from reverse transcription, PCR, and sequencing platform-specific errors (58). Together, these new complete genome sequences doubled the number available on GenBank.
Distinct lineages and multiple introductions of the EV-D68 viruses in Kansas City during the 2014 outbreak.Phylogenetic analyses of the EV-D68 VP1 gene and complete genome sequences revealed that the viruses that circulated in Kansas City during the 2014 outbreak did not form a single monophyletic group, such that EV-D68 was clearly introduced multiple times into the city (Fig. 2). More broadly, global VP1 phylogenies revealed the presence of the three major clades of EV-D68—A, B, and C—which have spread worldwide (Fig. 2A), consistent with previous studies (18, 38). Clades A and B could be further subdivided into two subclades each: A1 and A2 and B1 and B2. All clades and subclades were supported by bootstrap values over 80% in the ML approach and by posterior probabilities of 1.0 in the BMCMC approach.
Evolutionary history of the VP1 and complete genome sequences of EV-D68. (A and B) Global ML phylogenies estimated for the VP1 gene (A) and full genome sequences (B). Clades and subclades are indicated by colors and described in the trees. The bootstrap support values (percent) from 1,000 ML bootstrap replications and posterior probabilities from BMCMC analyses (10,000 tree samples) are shown for each clade/subclade with separation by a slash. Sequences reported by this study are shown in red and described in the key. Phylogenetic trees were rooted by the oldest EV-D68 sequence in GenBank, the Fermon strain collected in 1962 in California, USA. (C) Magnification of subclade B1 viruses from the VP1 phylogeny (A). Sequences generated in this study are marked by red squares in the tree. Sequences from acute flaccid myelitis (AFM) patients in the United States are marked by black squares. Sequences collected from outside the United States have taxon names in blue. The scale bars represent the numbers of nucleotide substitutions per site.
All 57 outbreak samples from Kansas City fell into subclade B1. Within this subclade, it was notable that the Kansas City viruses fell into two monophyletic groups and were interspersed among those viruses sampled from different states in the United States (New York, California, and Colorado) and from different countries (Canada or France), indicating frequent gene flow between populations (Fig. 2C). However, because of low bootstrap support values and frequent polytomies, it is difficult to define exact monophyletic groups of viruses from Kansas City on global VP1 phylogeny. Therefore, we estimated the number of independent viral introductions into the Kansas City population using a parsimony reconstruction of geographic transitions, utilizing 5,000 subsampled BMCMC topologies to reflect topological uncertainty. From this, we documented an average of 10 independent introductions of EV-D68 in to Kansas City during the outbreak, with a 95% confidence interval between 5 and 14 introductions.
The large and well-supported monophyletic group of viruses in subclade B1 (at the top of Fig. 2C) includes 2014 U.S. outbreak samples collected from New York, Missouri, and California; two 2014 samples from France; and one 2014 sample from Canada. Again, this phylogenetic pattern is indicative of both the national and international movement of viruses. Seven Kansas City sequences did not fall in this cluster; four were closely related to each other, while another sequence grouped with a virus from Canada. The remaining two Kansas City sequences formed separate branches in the tree, suggesting that they also represent independent introductions into the region. The remaining subclade B1 viruses are preoutbreak samples located near the root of the B1 cluster, including four U.S. isolates collected in 2013, a lone 2012 isolate from Italy, and, notably, three Chinese isolates sampled in 2011 (Fig. 2C). In addition, only five U.S. outbreak sequences collected in 2014 fell outside subclade B1: two from New York (NY73 and NY74), one from Kentucky (US.KY.14.18951), and one from Illinois (US.IL.14.18952), all of which fell in subclade B2. Sequence US.KY.14.18953 from Kentucky fell into subclade A2, while the preoutbreak sample collected in Nashville, Tennessee, in 2012 fell into subclade A1. Finally, there are two Fermon strain sequences in our analysis, one downloaded from GenBank and the prototype Fermon strain that we purchased from ATCC and sequenced. These two sequences clustered near the root, reflecting its sampling date in the 1960s (Fig. 2A).
The phylogeny of complete EV-D68 genome sequences presented the same picture of three major clades circulating worldwide and multiple viral introductions to Kansas City during the 2014 outbreak (Fig. 2B and 3). The nucleotide similarities among the clades ranged from 90.7% to 92.7%: clade A and clade B, 90.7%; clade A and clade C, 91.6%; and clade B and clade C, 92.7%. Nonstructural genes (2A, 2B, 2C, 3A, 3B, 3C, and 3D) exhibited levels of heterogeneity similar to those of structural genes (VP4, VP2, VP3, and VP1), with nucleotide similarity ranging from 96.5% to 97.1%. VP1 is the most variable gene across the genome as a whole, displaying 96.5% nucleotide similarity.
Phylogenetic trees of the complete genome and individual genes of EV-D68. Sequences in different clades are colored as described in the key. Phylogenetic trees inferred for different regions of the EV-D68 genome are indicated in black. Phylogenetic trees were rooted using the oldest EV-D68 sequence (the Fermon strain) available on GenBank, and scale bars represent the numbers of nucleotide substitutions per site.
EV-D68 recombination.We found no overt evidence of homologous recombination in the VP1 data set. However, a single recombinant virus—strain US.KY.14.18951—was identified in the complete genome data set, supported by P values of less than 1 × 10−15 in all methods. To locate the recombination breakpoints, a carefully chosen subset of sequence data (n = 22), including the recombinant strain and the closest parental lineages, was analyzed using similarity plots and Bootscan analysis in Simplot v3.5.1. This identified one significant recombination breakpoint in VP2, at nucleotide positions 1502 to 1576 (supported by a maximum χ2 sum of 34.3), that separated the genome into two nonrecombinant regions (Fig. 4) and which was strongly supported by phylogenetic trees inferred on either side of the breakpoint (Fig. 4C). Although the Bootscan analysis provided evidence for a second recombination event in the 3D gene region, this lacked statistical support. Hence, these data suggest that there was an intersubclade recombination event between B2 and B1. In contrast to subclade B1, which is the main U.S. 2014 outbreak lineage, the majority of subclade B2 sequences were collected from Asia and Europe during 2009 to 2014, with the exception of four U.S. 2014 outbreak sequences. Among the latter, one is the confirmed recombinant, and another is the potential parental strain, US.IL.14.18952 (Fig. 2A).
Recombination analyses of the full genome of EV-D68. (A) Bootscan analysis. A sliding window of 300 nt was utilized (step size of 10 nt), with neighbor-joining phylogenetic trees with 100 bootstrap replicates inferred in each case. Bootstrap values supporting the clustering of query sequence (i.e., US.KY.14.18951) with different groups of reference sequences (i.e., subclade B1, subclade B2, and others) were recorded. (B) Schematic diagram of the EV-D68 genome. (C) Phylogenies inferred for nonrecombinant fragments identified by Simplot and Bootscan (A). The putative recombinant strain is marked by the black square in the trees. Clades are indicated by the colors described in the key. Phylogenetic trees were midpoint rooted for clarity only.
No association between viral phylogeny and demographic and clinical metadata.We performed two phylogeny-trait association statistics (AI and PS) to identify whether EV-D68 outbreak strains might exhibit some phylogenetic clustering by disease severity. Notably, however, we found no significant association of age distribution, gender ratio, or clinical symptomatology with viral genetic variation (i.e., phylogenetic position) in the Kansas City cohort (Table 2). Hence, phylogenetic clustering by disease phenotype was not greater than that expected by chance alone. Similarly, although no EV-D68-associated AFM cases were reported in Kansas City, recently reported EV-D68 sequences collected from AFM patients from Colorado and California did not form a monophyletic group (38) but were interspersed with sequences of EV-D68 isolates from non-AFM patients (Fig. 2C; sequences from AFM patients marked by black rectangle).
Results of the phylogeny-trait association tests for particular demographic and clinical characteristics of EV-D68 patients within Kansas City, MO, 2014
DISCUSSION
We developed a high-throughput method to sequence complete genomes of EV-D68, from which we were able to obtain 59 from the major 2014 outbreak, including 57 from Kansas City. This doubles the number of publicly available genome sequences of EV-D68 and allowed us to identify both multiple introduction events into a single community and multiple cocirculating lineages worldwide, although there was no association between phylogenetic history and a range of demographic or clinical features. Notably, we also identified the first recombination event in EV-D68.
Until now, only limited full-length genome sequences of EV-D68 from the 2014 EV-D68 outbreak in the United States have been available in the public domain, with most coming from clinical isolates. With our high-throughput sequencing method, we were able to generate overlapping amplicons from primary specimens to obtain the complete viral genome with a 95% success rate. The success of this method compared to the single-amplicon approach is most likely due to RNA secondary structure in the 5′ region of the virus untranslated region that could impact the efficiency of RT-PCR.
Our phylogenetic analyses are notable in that they suggest a relatively complex molecular evolution of EV-D68 during the 2014 outbreak. Before this outbreak, EV-D68 was rarely reported and was associated with mild respiratory illness, although small outbreaks had been documented since 2005 (15, 16, 18–33). Previous studies showed that the A, B, and C clades circulated or cocirculated during different time periods in different geographic regions (15, 16, 18, 19, 21–32, 59). We further defined subclades A1, A2, B1, and B2. In line with our subclade definition, subclades A1 and B2 were endemic and found in many countries before the 2014 outbreak, including Thailand during 2005 to 2010 (26), the United Kingdom during 2009 to 2010 (27), China during 2009 to 2012 (24), Philippines from 2008 to 2014 (21–23), and the Netherlands from 1994 to 2013 (15, 29). Clade C was relatively geographically restricted, being reported in Japan during 2005 to 2010 and in Italy during 2010 to 2012 (16, 19, 30). During the 2014 outbreak, a novel subclade, B1, seems to have rapidly emerged and was dominant during the U.S. outbreak (38, 60, 61). During the same period, subclades A1 and B2 continued to be isolated in European countries (40, 41, 62–65). In addition, we found multiple viral introductions from other localities into the single locality of Kansas City during the 2014 U.S. outbreak despite the relatively small number of sequences collected (60), indicating that EV-D68 exhibits relatively fluid spatial dynamics within the United States.
Our study is noteworthy for the observation of recombination in EV-D68. Although recombination occurs frequently among members of the Enterovirus genus, with interspecies recombination documented between rhinovirus A and rhinovirus B and between EV-A and EV-B (33), and intraspecies recombination within both EV-A71 and EV-B81 (44, 66), it had not been previously definitely demonstrated to occur in EV-D68. Indeed, a previous report of phylogenetic incongruence between the 5′ UTR and VP1 could not be confirmed since only partial 5′ UTR and VP1 sequences were amplified (22). In contrast, our analysis of >100 complete EV-D68 genomes provided compelling evidence for recombination between outbreak subclade B1 and endemic subclade B2 strains, with a breakpoint in VP2. In recombination events in other enteroviruses, most breakpoints have been documented between the structural and nonstructural regions, or within the nonstructural region (44, 67). Interestingly, the previous report of phylogenetic incongruence between the 5′ UTR and VP1 in EV-D68 suggests the possibility of a recombination event within the structural region (22), although this clearly needs to be confirmed. The recombinant strain identified here, US.KY.14.18951 from Kentucky, and the potential parental strain in subclade B2, US.IL.14.18952 from Illinois, were both isolated in August 2014, at the start of the 2014 outbreak (68). In turn, the occurrence of recombination indicates that multiple distinct strains cocirculated at the beginning of the outbreak, as is also suggested by our phylogenetic analysis. The impact of VP2-centered recombination in EV-D68 evolution remains to be defined.
One of the most striking features of the 2014 outbreak was that most EV-D68 infections were associated with severe respiratory diseases in children (36, 40, 41, 60–65, 69). However, we found no association between age, gender, or clinical symptoms and different viral subclades in our Kansas City population. Hence, whether sequence data can predict the clinical outcome is still unclear, and our sample size was relatively small. One finding suggestive of an association between viral sequence diversity and clinical features is that AFM-associated strains from the United States during the 2014 outbreak were all members of the novel outbreak subclade B1 (38), although we were unable to identify viral genetic signatures of disease severity. Indeed, there are a number of features that argue against a strong strain basis to virulence. First, of the four cases of neurological conditions before the 2014 outbreak, two cases from Kenya (strain HEV126010 in 2010 and strain HEV199011 in 2011) were associated with viruses from subclade A1, the endemic subclade (31), while sequence information concerning the two U.S. patients is unavailable. Second, although all U.S. 2014 AFM cases fall in subclade B1 (38), two of three European AFM cases (40, 41) are from subclade B2. Third, the subclade B1-associated AFM sequences do not form a single cluster but rather are interspersed with those representing non-AFM cases (Fig. 2C). It was recently suggested that a G272 mutation in the 5′ UTR and a series of amino acid changes in the clade B1 polyprotein (T291 in VP2, A341 in VP3, N860 in VP1, N927 in 2A, K2005 in 3D, and G1108 at the 2B/2C junction) might have increased neurovirulence in EV-D68 (38). However, G272 was present in all the clades studied here and the polyprotein mutations were B1 specific, not AFM specific, as they were also found in non-AFM patients. In addition, two non-B1 AFM sequences, reported before the U.S. 2014 outbreak, possess none of the polyprotein mutations. Hence, the concentration of AFM cases in subclade B1 may simply reflect the predominance of B1 infections in the United States (i.e., a founder effect). Indeed, it is pertinent that EV-D68 viruses could not be isolated from cerebrospinal fluid of AFM patients from the 2014 U.S. outbreak (38).
Although we considered only a single city in the United States during the outbreak, these results are likely to be applicable to the spread of EV-D68 in other modern and highly mobile populations. For an improved understanding of the factors determining possible spatiotemporal differences in EV-D68 infection and transmission, continuous global monitoring of the clinical and molecular epidemiology of EV-D68 by representative surveillance systems should clearly be a public health priority.
ACKNOWLEDGMENTS
The clinical sample and data collection for this study were supported by a National Institute of Allergy and Infectious Diseases grant (AI U19-AI-095277). The sequencing work was supported by the NIAID/NIH Genomic Centers for Infectious Diseases (GCID) program (U19-AI-110819). Vanderbilt birth cohort and respiratory viral surveillance was supported by NIAID U19-AI-95227. E.C.H. was supported by an NHMRC Australia Fellowship (AF30).
The content is solely the responsibility of the authors and does not represent official views of the National Institutes of Health.
S.R.D., Y.T., R.S., E.C.H., T.V.H., and J.D.C. conceived this study. F.H., J.E.S., R.S., J.D.C., and T.V.H. collected clinical samples and clinical data and performed EV-D68-specific PCR. A.S., X.L., and R.A.H. performed RNA extraction, genome amplification, and viral sequencing. N.F. and T.B.S. assembled and analyzed the genomes and finished genome sequences as needed. Y.T., T.T.-Y.L., E.C.H., and S.R.D. analyzed the data. Y.T., T.T.-Y.L., E.C.H., and S.R.D. wrote the manuscript, and all authors reviewed and approved the final version.
We have no conflicts of interest to declare.
FOOTNOTES
- Received 18 September 2015.
- Accepted 30 November 2015.
- Accepted manuscript posted online 9 December 2015.
- Address correspondence to Suman R. Das, sdas{at}jcvi.org.
Citation Tan Y, Hassan F, Schuster JE, Simenauer A, Selvarangan R, Halpin RA, Lin X, Fedorova N, Stockwell TB, Lam TT-Y, Chappell JD, Hartert TV, Holmes EC, Das SR. 2016. Molecular evolution and intraclade recombination of enterovirus D68 during the 2014 outbreak in the United States. J Virol 90:1997–2007. doi:10.1128/JVI.02418-15.
REFERENCES
- Copyright © 2016, American Society for Microbiology. All Rights Reserved.
