Intra- and interpatient evolution of enterovirus D68 analyzed by whole-genome deep sequencing

Abstract

Worldwide outbreaks of enterovirus D68 (EV-D68) in 2014 and 2016 have caused serious respiratory and neurological disease. To investigate diversity, spread, and evolution of EV-D68 we performed near full-length deep sequencing in fifty-four samples obtained in Sweden during the 2014 and 2016 outbreaks. In most samples, intrapatient variability was low and dominated by rare synonymous variants, but three patients showed evidence of dual infections with distinct EV-D68 variants from the same subclade. Interpatient evolution showed a very strong temporal signal, with an evolutionary rate of 0.0039 ± 0.0001 substitutions per site and year. Phylogenetic trees reconstructed from the sequences suggest that EV-D68 was introduced into Stockholm several times during the 2016 outbreak. Putative neutralization targets in the BC and DE loops of the VP1 protein were slightly more diverse within-host and tended to undergo more frequent substitution than other genomic regions. However, evolution in these loops did not appear to have been driven the emergence of the 2016 B3-subclade directly from the 2014 B1-subclade. Instead, the most recent ancestor of both clades was dated to 2009. The study provides a comprehensive description of the intra- and interpatient evolution of EV-D68, including the first report of intrapatient diversity and dual infections. The new data along with publicly available EV-D68 sequences are included in an interactive phylodynamic analysis on nextstrain.org/enterovirus/d68 to facilitate timely EV-D68 tracking in the future.

Keywords: enterovirus; interpatient evolution; intrahost evolution; whole-genome deep sequencing.

Figure 1.

Example of coverage across the EV-D68 genome in a representative clinical sample. Note that the coverage varies somewhat in different parts of the genome due slightly different quantity of amplified DNA from the four overlapping amplicons in the material that was used for sequencing library preparation. The four overlapping are indicated by the horizontal gray bars.

Figure 2.

Accuracy of minor variant frequencies. Panel (A) shows the consistency of minor variant frequencies across two replicate extraction, RT-PCR, and sequencing of two samples. Minor variation was relatively consistently recovered down to a frequency of about 1 per cent. Panel (B) shows the distribution of the frequencies of non-consensus calls across all sites covered in excess of 2,000-fold. At the majority of sites, no variation above 0.001 is observed. Note that this variation includes within sample variation, RT-PCR errors, and sequencing errors.

Figure 3.

Within sample diversity. Panel (A) shows the distribution of the number of variable sites in coding regions among different coding positions. At very low-frequency cutoffs, variable sites are approximately equally distributed in codons. The distribution rapidly changes once the frequency cutoff exceeds 0.3 per cent, beyond which variable sites are found predominantly at third codon positions and more rarely at first and second codon positions. Panel (B) shows inverse cumulative distribution of the number of samples that have more than a certain number of sites that are variable above a level of 3 per cent. Separate distributions are shown for all sites and for first, second, and third codon position in the coding region of the genome.

Figure 4.

EV-D68 intrapatient variability across the genome and codon positions. The three rows show how intrapatient variability is distributed across the EV-D68 genome for first, second, and third codon position, respectively. The panel on the left show iSNVs along the entire polypeptide, the panels on the right zoom into VP1 and highlight the BC and DE loops. This figure includes data from one sample per patient (the first sample, for patients sampled twice), but excludes the three dual infected samples (see main text).

Figure 5.

Minor iSNVs across the EV-D68 genome for three samples with putative dual infections. These three samples had more than twenty highly covered sites with minor variants in excess of 3 per cent frequency. On the left, iSNVs are colored by codon position or as non-coding, and given opacity depending on their frequency. The majority of these variants were at third positions and had similar frequencies across the amplicons, indicated by gray bars in the lower left panel. iSNVs in three amplicons (Amplicon 1 of sample SWE_046 and Amplicons 2 and 3 of sample SWE_007), however, were found at substantially lower frequency, possibly due to primer mismatches. The three panels on the right are indexed by iSNV order on the genome and show linkage disequilibrium between iSNVs (>1%) close enough to each other that they were covered at least 100-fold by the same sequencing read. Almost all of these variants at in complete linkage (dark red). A number of iSNVs just above 1 per cent in sample SWE_045 and SWE_046 are likely variants in the background of the dominant variant. Those are in complete ‘anti-linkage’ with neighboring iSNVs are ∼10 per cent (dark blue).

Figure 6.

Phylogenetic analysis. The left panel shows a phylogenetic tree as rendered by nextstrain with the sequences labeled by region. Those from Europe, including those from Sweden, are colored in orange. The right panel is a zoomed-in view of the boxed area on the left-hand panel colored by country and highlighting the Swedish sequences from the 2016 outbreak (orange). Swedish sequences often cluster together, but are interspersed with sequences from Canada (dark red) and the USA (red), implying multiple introductions that then spread locally.

Figure 7.

Temporal signal. A scatter plot of root-to-tip distance vs. time indicates an extremely strong temporal signal with an evolutionary rate of μ = 0.0039 ± 0.0001 substitutions per site and year.

Figure 8.

Interpatient diversity of EV-D68. Amino acid substitution events on the phylogeny in Fig. 6 left were enumerated using nextstrain. Top panel shows the number of such events for each codon along the entire coding sequence. Bottom panel shows events in VP1, where peaks were observed in the BC and DE loops, as well as both ends of the protein.