Enterovirus D68 (EV-D68) has been increasingly recognized as an emerging virus, with outbreaks occurring primarily in children (Piralla et al., 2018; Kamau et al., 2019; Andrés et al., 2022). First reported in 1962 in California, EV-D68 was sporadically reported worldwide until 2013, with small outbreaks described in the Philippines (2009–2011; Imamura and Oshitani, 2015), Japan (2010; Hasegawa et al., 2011), United States of America (USA; 1970–2005; Khetsuriani et al., 2006) and in the Netherlands (2010; Rahamat-Langendoen et al., 2011). The largest known outbreak to date was reported in the USA in 2014, with over a thousand cases in children presenting with severe respiratory disease (Midgley et al., 2015). During the 2014 outbreak, a number of EV-D68 cases were linked to acute flaccid myelitis (AFM), presenting with asymmetric flaccid limb weakness and cranial nerve dysfunction (Messacar et al., 2015; Vogt and Crowe, 2018). In the same year, EV-D68 was also reported in Europe, Southeast Asia and Chile (Holm-Hansen et al., 2016). Since 2014, EV-D68 outbreaks have occurred in a biennial pattern between 2016 and 2018, typically in even years in temperate climates (Howson-Wells et al., 2022) and emerging most frequently in the USA, Europe, Argentina and Taiwan (Knoester et al., 2017; Hu and Chang, 2020). While EV-D68 typically peaks in late summer and autumn in temperate climates (Hodcroft et al., 2022), an off-season upsurge was observed across Europe during the 2019–2020 winter season (Midgley et al., 2020).

EV-D68 has been shown to have substantial genetic diversity (Kramer et al., 2018; Howson-Wells et al., 2022). Sequences can be separated into three main clades: A, B, and C, which can be further divided into A1–A2 and B1–B3 subclades (Hodcroft et al., 2022). The EV-D68 genome is approximately 7,500 nucleotides in length and encodes a single polyprotein which consists of four structural proteins (viral proteins 1–4) and seven non-structural proteins (2A–C and 3A–D; Eshaghi et al., 2017; Dyrdak et al., 2019). Generally, EV-D68 sequences have been investigated using the hypervariable viral protein 1 (VP1) gene, which is present on the viral capsid and has been used as a target for subtype differentiation (Harvala et al., 2018). Since 2014, subclades B1 and B3 have become dominant (Hodcroft et al., 2022). The emergence of the A2 subclade was additionally reported in East Asia in 2016 (Midgley et al., 2020). By 2018, clusters of subclades B3 and A2 appear to be predominating worldwide.

Whole-genome sequencing (WGS) has been used previously to investigate the diversity and evolution of EV-D68 during the 2014 and 2016 outbreaks (Dyrdak et al., 2019). Different subclades have been identified in respiratory and neurological samples, with varying rates of amino acid substitutions observed (Zhang et al., 2016; Wang et al., 2017; Hodcroft et al., 2022). In this study, we obtained near-complete EV-D68 sequences from patients presenting with respiratory or neurological disease at a regional university hospital in the Netherlands during the winter season of 2019–2020. Using WGS, we aimed to explore our patient cohort, along with the circulating clades responsible for the rise in cases with the intention of contributing to a greater understanding of potential evolutionarily changes within EV-D68.

In this study, we provide clinical and phylogenetic information from hospitalized patients from an off-season EV-D68 upsurge in the winter season of 2019–2020. Using a modified WGS approach, we were able to obtain near-complete EV-D68 genomes from a wide range of Ct values, including Ct values >30 and sample material, including CSF. The summer/autumn seasonal circulation is well recognized for all enteroviruses in temperate climates, with a biennial recurrence of EV-D68 observed in North America and Europe (Elrick et al., 2021; Fall et al., 2022). Interestingly, our data shows a general shift towards the later months, with the upsurge emerging in an odd year (Figure 1). This off-season occurrence was also observed by other European countries (Midgley et al., 2020). Although the seasonality of enteroviruses is likely determined by environmental factors such as temperature and air humidity (Fares, 2013; Pons-Salort et al., 2018), the biennial occurrence of EV-D68 is thus far not understood. While factors including herd immunity, human behavior and circulation of other (entero-) viruses may play a part, the role and occurrence of viral variants requires more investigation.

The majority of patients in our study had underlying conditions, with prematurity highest among children and typically multiple comorbidities in adults, echoing previous studies (Lau et al., 2016; Fall et al., 2022). Overall, we observed 13 cases (65%) of respiratory illness and three cases (15%) of neurological illness (all in children; Table 1), two of which had AFM and one had EV-D68 meningitis. In our small cohort, we found a relatively high number of cases requiring admission to hospital, with 75% of patients admitted for over 2 days. This most likely reflects the tertiary-care character of our hospital. A similar ratio of severe illness was additionally observed within Europe during the same time period (Midgley et al., 2020). Interestingly, we identified one A2 subclade from a child (<1 year) with meningitis. Although detection of EV-D68 in CSF is rare (Fares, 2013), to the best of our knowledge, we report the first detection of an A2 subclade associated meningitis in CSF. EV-D68 was also detected from one fecal sample from a patient presenting with respiratory and diarrhea illness. Although EV-D68 detection from fecal samples is uncommon, it could provide an alternative for widespread population surveillance studies (Corpuz et al., 2020; Tedcastle et al., 2022). It is important to note that the true burden and circulation within the community is difficult to estimate, given that community surveillance is not performed and subtype differentiation is typically only performed for severer cases following clinical consultation.

Genotyping enteroviruses is important for understanding changes in epidemiology, which could aid outbreak preparedness and patient management (Harvala et al., 2018). WGS has been previously used to explore viral dynamics, particularly during outbreak scenarios and especially during the SARS-CoV-2 pandemic (Oude Munnink et al., 2021). We adapted a previously applied EV-D68 WGS method for ONT platforms to characterize our upsurge samples, including previously untypeable samples.

Using the Nextstrain pipeline, we could perform more real-time tracking of EV-D68 evolution and spread overtime (Figure 1). A high genetic diversity can be observed among the available EV-D68 sequences, indicating diversification and persistence even in low circulation years, similarly reported by Hodcroft et al. (2022). Interestingly, there appears to be continuous emergence in conjuncture with clade replacement prior to 2018, with the rise of the predominating B3 subclade. In our study, 18 samples clustered with the B3 subclade within three different branches, suggesting introduction at several different timepoints in the Netherlands, followed by local spread. Subclade B3 exhibits high levels of geographical circulation (and mixing) since its emergence in 2016 (Wang et al., 2017), having been found in several countries globally (Midgley et al., 2020). The detection of the A2 subclade additionally suggests an ongoing low-level (and potentially undetected) transmission in Europe. Given their increasing observation between 2014 and 2020 (Lau et al., 2016; Midgley et al., 2020; Hodcroft et al., 2022), it could render this finding potentially important for monitoring. Moreover, although the majority of sequences from our study appear to cluster with USA strains, as some countries deposit more sequences than others, it could also result in a geographical bias. Our time-scaled tree generated an estimated evolution rate of 3.05 × 10⁻³ for the near-complete genome, which is slightly lower than 3.8 × 10⁻³ in 2018 (Dyrdak et al., 2019), but higher than another estimate of 2.99 × 10⁻³ from 2017 (Eshaghi et al., 2017). Some patient sequences (patient 1, 4,12) appeared to have higher genetic diversity than average for the sampling date (Supplementary Figure S3; Supplementary Table S7), suggesting a possible higher level of viral evolution within these patients

The generation of high-quality near-complete genomes potentially allows the investigation of inter-and intra-host relations of SNV (Dyrdak et al., 2019; Cassidy et al., 2022). As seen in the SARS-CoV-2 pandemic, variants may explain changes in transmissibility and pathogenicity (Yang et al., 2022). We report a high level of variation within our study, particularly in the VP1 and 3C (transcribing the EV-D68-encoded protease) genes, however in this small cohort, any link between SNV and clinical outcome would not be statistically relevant. Moreover, coverage bias may have been introduced in our study with some samples generating a low volume of reads, while other samples obtaining sequencing depths higher than 100,000× in some regions. This coverage bias was particularly prominent in regions with overlapping primer fragments (Supplementary Figure S2). As a result, we applied a coverage cut-off of 100,000× to reduce the chances of sequencing or PCR errors which could be called wrongly due to the sheer number of reads. Additionally, although continuous improvements in ONT error rate have been made over the recent years (Wick et al., 2021), the overall high error rate compared to Illumina sequencing could also be a limitation. Variants with the highest frequency (in 90% of patient sequences) were subsequently revealed to be highly conserved within the B3 and A2 subclades (>90% in our reference database).

The high number of variants within the hypervariable VP1 gene could highlight the impact of immune pressure and its continual influence to shape this region (Opanda et al., 2014; Eshaghi et al., 2017). Other studies have identified rapidly evolving neutralizing epitopes in the VP1 gene, particularly in the BC and DE loops, with higher amino acid substitution rates than the rest of the genome (Hodcroft et al., 2022). With the emergence of new EV-D68 strains, given the high genetic diversity observed over time (Figures 2, 3), coupled with declining immunity (Kramer et al., 2018), it could indicate a growing susceptible population. At the same time, the requirement of re-amplification in some sample fragments could increase the chances for PCR driven mutations (Salk et al., 2018).

Finally, we show evidence of high levels of episodic positive diversifying selection events within the EV-D68 population and within our patient dataset scattered along the genome, particularly in the 3C and VP1 genes (referred to as “hotspot” regions; Muslin et al., 2019). As our patient samples are the most recently deposited (near-full) EV-D68 sequences, it could highlight the areas which are most likely being selected in the future. Interestingly, the 3C gene has a central role in inhibiting antiviral immunity (Xiang et al., 2015). It has been speculated that the now dominant B3 subclade has arose from diversifying selection and likely represents a newly emerging strain, as opposed to evolving from the B1 or B2 clades (Wang et al., 2017; Hodcroft et al., 2022). Given the strong selection pressure observed (in this case positive diversifying selection) of subclade B3, along with a potentially increasing naïve population, it could have driven its predominance since 2018.

We report an off-season EV-D68 upsurge likely caused by highly diverse subclades and associated with severe respiratory and neurological disease, particularly in children and immunocompromised adults. It could be that the positive episodic or diversifying selection, along with persistent yet undetected circulation may be driving EV-D68 evolution, however at this point this is still a working hypothesis. Furthermore, high levels of non-synonymous mutations, particularly in the surface proteins, could suggest monitoring closely for potential challenges with routine Sanger sequencing. Applying WGS revealed an dynamic picture of EV-D68, highlighting evolutionary patterns. With the re-emergence of EV-D68 following the COVID-19 lockdown and poliovirus recently detected in non-epidemic countries, it is important to closely monitor this virus which may become a threat to public health.

The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA900246.

HC: study design, laboratory and data analysis, writing-original draft preparation, and writing-reviewing and editing. EL-F: data analysis and writing-reviewing and editing. LS: study design, laboratory and data analysis, and writing-reviewing and editing. CL-B and HGMN: conceptualization, editing, and supervision. All authors contributed to the article and approved the submitted version.

Hayley Cassidy and Leonard Schuele received funding from the European Union’s Horizon 2020 research and innovation program, under the Marie Sklodowska-Curie grant agreement 713660 (MSCA-COFUND-2015-DP “Pronkjewail”). Erley Lizarazo-Forero was funded by Quality Control Molecular Diagnostics (QCMD, Glasgow, Scotland) under an unrestricted grant.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.