Data Collection

An invitation to participate in this study was sent on November 13, 2018, to co-authors of the E30–2018 study (2) through the European Centre for Disease Prevention and Control (ECDC) Epidemic Intelligence Information System Vaccine-Preventable Diseases platform (https://www.ecdc.europa.eu/en/publications-data/epidemic-intelligence-information-system-epis), and to members of the European Non-Polio Enterovirus Network (ENPEN; https://www.escv.eu/enpen). We requested pseudonymized data from 2016–2018 with sample identifier, sampling date, specimen type, and sequence in FASTA be sent to ECDC secure file transfer protocol server by January 7, 2019. We also collected optional data, such as patient age, clinical presentation, whether they were hospitalized, and infection outcome. We excluded submissions without virus sequence data (Appendix Figure 1).

Sequence Data Collection

We requested that the FASTA sequence data contain the VP1 gene and collected 1,784 records (Appendix Figure 1). Sequences were obtained from enterovirus-positive samples by using 5′ UTR PCR (23) and typed within the VP1 gene by using Sanger sequencing (2,24). We excluded sequences with indicators of poor sequence quality, such as >2 ambiguous or undefined bases, in-frame stop codons, identical to reference E30 strains; sequences of the wrong type, such as E3; or sequences shorter than 200 basepairs or spanning a non-VP1 region. In total, we had 1,407 study sequences that comprised 2 nonoverlapping regions, 1,262 sequences from region 1 (nt positions 2543–2902, according to the prototype E30 strain Bastianni, GenBank accession no. AF311938) and 145 sequences from region 2 (nt positions 2916–3428). Of these, 1,329 sequences were collected during 2016–2018 and 78 during 2010–2015. We used the 2010–2015 sequences for phylogenetic reconstruction but excluded these from further data analysis (Appendix Figure 1).

For additional analysis of the 3D polymerase (3Dpol) region, we randomly selected records from each clade to ensure fair distribution of sequence data. We asked participants to send either extracted RNA in a QIAGEN (https://www.qiagen.com) spin column at room temperature for next-generation sequencing (NGS) or to conduct 3Dpol sequencing of the 549 nucleotides, as previously described (17).

Epidemiologic and Statistical Analyses

We descriptively analyzed clinical symptoms and age. Patients were stratified into the following age groups: <3 months, 3–23 months, 2–5 years, 6–15 years, 16–25 years, 26–45 years, 46–65 years, and >65 years. Crude odds ratios with 95% CI were used to express magnitude of association between continuous or categorical variables in multivariate logistic regression.

Next-Generation Sequencing

Stool suspensions and CSF samples were processed to remove as much nonviral material as possible by using centrifugation, filtration, and endonuclease treatment. RNA was extracted by using the MagNAPure 96 (Roche Diagnostics, https://www.roche.com) automated extraction kit or QIAGEN filters and eluted in 50 µL of elution buffer (Appendix).

Complementary DNA (cDNA) and double stranded DNA (dsDNA) were generated and purified (Appendix). For tagmentation and library preparation, the Nextera XT DNA Library Preparation Kit (Ilumina, https://www.illumina.com) was used according to the manufacturer’s instructions. Runs were performed on the Nextseq (Ilumina). Raw data were processed by using Jovian (D. Schmitz et al., unpub. data, https://github.com/DennisSchmitz/Jovian) (Appendix).

Nucleotide Sequences and Phylogenetic Analysis

We conducted VP1 phylogenetic reconstruction with the 1,407 study sequences and 324 sequences extracted from Genbank. We selected region 1 for clade analysis because it is more commonly used for enterovirus typing (24). We performed analysis of region 2 sequences based on sequence clustering, in which both region 1 and 2 were available, such as full-length sequences or sequences spanning the entire VP1 gene. Sequencing of the 540 nt 3Dpol gene, positions 5825–6364, also was provided for 12 samples with region 1 sequences (Appendix Figure 1). Sanger sequencing indicated that samples did not display double infection and that VP1 and 3Dpol were from 1 virus. Complete genomes (≈7.3 kb) were generated for 48 sequences by using NGS. To compare 3Dpol groupings within EV species B, we downloaded all sequences available from GenBank as of October 18, 2019, that were >70% complete between positions 5825–6364 with <6 ambiguous base positions and <6 undetermined bases and without stop codons. We aligned the downloaded sequences with complete genomes or 3Dpol sequences from our study.

We aligned data by using sequence editor SSE version 1.3 (http://www.virus-evolution.org). We generated maximum-likelihood and neighbor-joining trees for VP1 and 3Dpol regions by using MEGA version 7 (https://www.megasoftware.net) with the optimal model (general time reversible plus invariant sites plus gamma distribution for rates over sites) and 100 bootstraps (25). We analyzed the species B dataset with neighbor-joining and maximum composite likelihood distances.

Nextstrain VP1 Phylodynamic Analysis

The dataset used for Nextstrain phylodynamic analysis comprised 1,285 sequences; 1,215 study sequences (region 1) and 70 complete VP1 sequences extracted from GenBank (Appendix Figure 1). We excluded sequences shorter than 250 bp, sequences from samples collected before 1958, and sequences deemed as outliers during phylogenetic reconstruction. We deemed these outliers recombinants with possible recombination breakpoints within the sequence fragment used made phylogenetic reconstruction impossible.

We aligned sequences by using MAFFT (26). We inferred a phylogenetic tree by using IQ-TREE (27) and generated time-resolved trees by using TreeTime (28) by estimating the mutation rate. When available, we attached to sequences data on country, sample type, E30 clade, age groups, and clinical data, such as whether patients were hospitalized and their symptoms. We provided the resulting Nextstrain build for viewing (https://nextstrain.org/community/enterovirus-phylo/echo30-2019/vp1).

Nextstrain VP1:3Dpol Tanglegram Phylodynamics

We used a dataset of 110 sequences to conduct 3Dpol analysis, including 48 complete genome sequences and 12 3Dpol sequences generated in this study and 50 sequences extracted from GenBank (Appendix Figure 1). We aligned 3Dpol sequences to the E30 reference sequence (GenBank accession no. MK238483) and inferred phylogenetic and time-resolved trees as we did for VP1, but we used a fixed clock rate of 4 × 10⁻³ substitutions/site/year during the time-resolved tree reconstruction. We provided the resulting Nextstrain tanglegram build for viewing (https://nextstrain.org/community/enterovirus-phylo/echo30-2019/3D:community/enterovirus-phylo/echo30-2019/vp1) and the codes for both VP1 and 3Dpol analyses (https://github.com/enterovirus-phylo/echo30-2019).

Genbank and ENA Accession Numbers

We deposited VP1 and complete genome sequences in GenBank under accession nos. KC309427–37, KY986976–7033, MK251835–6, MK372854–80, MK507733–7, MK814991–6288, and MK895104–9 and 3Dpol sequences under accession nos. MN395293–303. We deposited NGS fastq reads in European Nucleotide Archive database under accession nos. SAM17101211–58.

Molecular Epidemiology and Demographics

Figure 1

Figure 1. Geographic distribution of echovirus 30 (EV30) clades, Europe, 2016–2018. Clades G1–G6 were detected among 1,329 EV30 cases from 22 countries. A) 2016; B) 2017; C) 2018.

Figure 2

Figure 2. Monthly distribution of echovirus 30 (EV30) clades G1–G6 detected among 1,329 sequences submitted from 22 countries in Europe during 2016–2018.

During 2016–2018, a total of 1,329 E30 records representing 1,292 cases that fulfilled the study criteria were submitted from 22 countries (Table 1; Appendix Figure 1). During those 3 years, the total number of E30 cases steadily increased (Table 1). The numbers varied per country per year, and we noted a clear upsurge in 2018 in several, but not all, countries (Table 1; Figure 1). Of the 1,329 records analyzed, 443 (33%) were from the United Kingdom; the Netherlands submitted 198 (15%) and Spain 162 (12%) records. Other countries submitted from 1 (<1%) to 117 (9%) records. Specimen type was reported for 1,312 (98.7%) records. Most (70%; 924/1,312) samples were cerebrospinal fluid specimens, but other specimen types included 269 (21%) from feces specimens, 102 (8%) from respiratory, and 17 (1%) from blood. During the study period, E30 records were submitted more frequently in summer months; 18.4% (n = 244) were submitted in June, 17.6% in July (n = 234), and 11.7% in August (n = 155) (Figure 2).

Age was available for 1,080 (83.6%) cases and ranged from 0 to 73 years with a mean age of 18.7 years. Children <5 years of age (n = 360, 33.3%) and adults 26–45 years of age (n = 409, 37.9%) were most affected (Table 2). Infants <3 months of age also were heavily affected (n = 223 cases, 20.6%) (Table 2).

Clinical information was available for 734 (56.8%) E30 cases, of which 380 cases had unknown symptomology. For most (28.7%, n = 211) cases, the recorded signs and symptoms suggested meningitis. Symptoms of acute flaccid paralysis were reported in 1 case, encephalitis in 3 cases, and meningoencephalitis in 8 cases. Fever, either as sole symptom or in combination with other signs and symptoms, was recorded in only 52 (7.1%) cases. Unfortunately, not all records were filled in completely, and clinical data were absent for some samples. Other signs and symptoms mentioned were gastrointestinal symptoms in 6 cases, respiratory symptoms in 6, rash in 2, other neurologic symptoms in 4, or other unspecified in 53 cases; 8 cases had no symptoms. We created an interactive representation of age and clinical features of sequences from E30 cases, which we made available on Nextstrain (https://nextstrain.org/community/enterovirus-phylo/echo30-2019/vp1).

Hospitalization status was available for only 17.6% (n = 228) of cases, only 5 of which had no hospitalization. The low fraction of hospitalization reported limited further analysis. No deaths were reported.

E30 Phylodynamics

Figure 3

Figure 3. Phylodynamic analysis of region 1 in a curated study of echovirus 30 (E30) viral protein 1 (VP1) sequences from 22 countries in Europe, 2010–2018. We constructed the bootstrapped maximum likelihood...

Among the 1,329 curated VP1 study sequences, 1,019 (76.7%) could be subdivided into 5 distinct clades, G1–G5, that showed >5% sequence divergence from one another (Figure 3, panel A). The mean divergence between VP1 nucleotide sequences of G1–G5 was 12.4%–15.2%, which translated to 2.8%–3.9% amino acid sequence divergence. Most (704, 53%) sequences belonged to G1, but 229 (17.2%) were in G2, 59 (4.4%) in G4, 4 (0.3%) in G3, and 2 (0.2%) in G5. These sequences all were assigned to GGII, 1 of 2 previously reported genogroups (7). The remaining 310 VP1 sequences formed a single clade, G6 (Figure 3), showing 20.6% mean nucleotide differences and 8.5% amino acid differences from the VP1 sequences within G1–G5 clades. G6 was sufficiently divergent from G1–G5 (GGII). The divergence falls within the nucleotide divergence between GGI–GGII (19%–22%) (7), and G6 can be considered a third genogroup, GGIII.

The phylogeny showed a rapid turnover of E30 clades over the 3 years sampled, shifting from G2 dominating in 2016 to G1 and G6 dominating in 2017 and 2018 (Figure 1). In 2016, 58.5% of strains were G2, and this genotype was identified in 11/22 (50%) countries. G2 was detected in only 8 countries in 2017 and only 4 countries in 2018. Similarly, G4 disappeared during 2016–2018. In 2016, both G1 and G6 were detected, G1 in 24.6% (n = 80) of virus strains in 10 countries and G6 in 6.5% (n = 21) of virus strains in 7 countries. Rates of detection for G1 and G6 steadily increased in 2017; G1 was detected in 59.2% (n = 292) of virus strains in 12 countries and G6 in 26.4% (n = 130) of virus strains in 7 countries (Figure 1). During the 2018 upsurge, 64.6% (330) of sequences reported in 17 countries belonged to G1, and 33.9% (173) in 11 countries belonged to G6 (Figure 1). These data indicate the occurrence of >2 distinct viruses dominating the upsurge in 2018 (Figure 1).

We used Nextstrain to create an interactive phylodynamic tree and map to explore relationships of the E30 study VP1 sequences in G1–G6 (https://nextstrain.org/community/enterovirus-phylo/echo30-2019/vp1) (Figure 3, panel B; Appendix Figure 2). We deemed G5 sequences as outliers and did not include these during phylogenetic reconstruction. Molecular clock analysis of the VP1 region revealed an estimated substitution rate of 5.12 × 10⁻³ substitutions/per site/per year, comparable to rates previously determined for a range of enteroviruses.

Most E30 G1 viruses were detected among infants <3 months of age (135/568, 24%) and in young adults 26–45 years of age (227/568, 40%) (Table 2). G2 (100/145, 68%) and G4 (29/59, 52%) were most frequent among children 3 months–15 years of age. G6 mainly was detected among children 3 months–15 years of age (100/308, 32.5%) and in adults 26–45 years of age (134/308, 43.5%). Only 2 cases of G3 and 1 of G5 were reported with age information (Table 2).

Amino Acid Diversity

Most E30 VP1 sequences within clades G1, G2, G4, and G6 displayed specific amino acid substitutions. G6 sequences predominantly displayed amino acid changes at position 56 (Y-F), position 84 within the BC loop (V-A), position 87 within the BC loop (E-D), and position 145 (V/I) compared with G1, G2, and G4. Most G1 and G6 sequences had a valine at positions 54 and 120 compared with the G2 and G4 sequences, which had an isoleucine. At position 122, most G4 sequences contained a leucine, whereas G1, G2, and G6 sequences contained a phenylalanine. Interactive data are available on Nextstrain (https://nextstrain.org/community/enterovirus-phylo/echo30-2019/vp1).

Recombination Analysis

Figure 4

Figure 4. Neighbor-joining tree of 3D polymerase (3Dpol) sequences of echovirus 30 (E30) study samples and sequences from previously described E30 strains. The tree was constructed from Jukes-Cantor corrected nucleotide sequence distances...

Figure 5

Figure 5. Tanglegram of echovirus 30 (E30) phylogenetic virus protein 1 (VP1) (right) and 3D polymerase (3Dpol) (left) by year of sample collection. We used 110 sequences and rendered the tanglegram by...

We used 110 sequences containing both VP1 and 3Dpol region and complete genome sequences to analyze recombination events between VP1 and the 3′ distal end of the E30 genome (Appendix Figure 1). The E30 3Dpol sequences formed a series of separate clusters interspersed with those of other species B types, indicative of many within-species recombination events during their diversification (Figure 4). We took the entire published sequence dataset and used a nucleotide sequence distance threshold of 8%, based on pairwise sequence comparisons, which divided sequences into distinct groups (Appendix Figure 3), comparable to those derived from a previous analysis of E30 RFs (17). Accordingly, species B could be divided into ≈442 RFs, an indication of the frequency and complexity of recombination events occurring during the evolution of this species. We used Nextstrain to generate an interactive tanglegram of VP1 and 3Dpol RFs (https://nextstrain.org/community/enterovirus-phylo/echo30-2019/3D:community/enterovirus-phylo/echo30-2019/vp1) (Figure 5).

We found that 3Dpol sequences of G1–G6 formed 8 recombination groups, which were separated by other published E30 variants and by other species B types (Figure 5). We noted that viruses within most VP1 clades were monophyletic in 3Dpol, but that G1 and G4 each had undergone further recombination (Figure 5; Appendix Figure 4), a split corresponding to the sublineages evident in the VP1-based tree. The split was identified as a time-related phenomenon, with G1 circulating in 2018 representing a different RF from G1 circulating during 2016 and 2017.